Primero, cargamos las librerías.

library(ggplot2)
library(ggpubr)
library(summarytools)
library(psych)
library(dplyr)
library(purrr)
library(readxl)
library(scatterplot3d)
library(reshape2)
library(MVN)
library(caret)
library(MASS)
library(biotools)
library(klaR)
library(cluster)
library(ggdendro)
library(gridExtra)
library(tidyverse)
library(factoextra)
library(gridExtra)
library(corrplot)

Importamos los datos. Como nuestro excel tiene tres pestañas, lo separamos en 3 data frames distintos. Datos1 para los datos del estudio 1, Datos2 para los del estudio 2 y Datos3 para los del estudio 3.

Datos1<- read_excel("C:\\Users\\Laura\\Desktop\\Datos arreglados.xlsx",1)
#View(Datos1)
Datos2<- read_excel("C:\\Users\\Laura\\Desktop\\Datos arreglados.xlsx",2)
#View(Datos2)
Datos3<- read_excel("C:\\Users\\Laura\\Desktop\\Datos arreglados.xlsx",3)
#View(Datos3)

Análisis exploratorio univariante

Datos ausentes (NA). Primero vemos qué datos son los que faltan. Observando las bases de datos, nos damos cuenta de que en el segundo estudio faltan BMC, BMD y longitud de las ratas con ID 222 y 223, y Tb.Sp, Tb.N, BMC,BMD, longitud, GPSurface, GPVolume de las ratas con ID 295 y 298. Tenemos que decidir si quitamos dichas ratas o simplemente omitimos los datos ausentes. Por ahora voy a trabajar omitiendo los datos ausentes. Vamos a realizar un análisis descriptivo separado por estudios.

ESTUDIO 1

Peso inicial

descr(Datos1[,14])
## Descriptive Statistics  
## Datos1$PesoInicial  
## N: 56  
## 
##                     PesoInicial
## ----------------- -------------
##              Mean        120.13
##           Std.Dev          8.75
##               Min        102.30
##                Q1        113.25
##            Median        118.95
##                Q3        127.35
##               Max        140.80
##               MAD         10.30
##               IQR         14.05
##                CV          0.07
##          Skewness          0.22
##       SE.Skewness          0.32
##          Kurtosis         -0.86
##           N.Valid         56.00
##         Pct.Valid        100.00
pesosiniciales1<-Datos1[[14]]
par(mfrow = c(1, 2),cex=0.7)
hist(pesosiniciales1, col="pink",xlab="Pesos iniciales", 
     ylab="Frecuencia", main="Distribución de los pesos iniciales estudio 1")
plot(density(pesosiniciales1), main="Densidad pesos iniciales estudio 1",
     xlab="Pesos iniciales", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(pesosiniciales1, main="Pesos iniciales estudio 1", outline=TRUE)

Peso intermedio

descr(Datos1[,15])
## Descriptive Statistics  
## Datos1$PesoIntermedio  
## N: 56  
## 
##                     PesoIntermedio
## ----------------- ----------------
##              Mean           184.64
##           Std.Dev            49.45
##               Min           139.20
##                Q1           158.65
##            Median           167.50
##                Q3           181.00
##               Max           334.00
##               MAD            16.53
##               IQR            22.08
##                CV             0.27
##          Skewness             1.83
##       SE.Skewness             0.32
##          Kurtosis             2.01
##           N.Valid            56.00
##         Pct.Valid           100.00
pesosintermedios1<-Datos1[[15]]
par(mfrow = c(1, 2),cex=0.7)
hist(pesosintermedios1, col="pink", xlab="Pesos intermedios", 
     ylab="Frecuencia", main="Distribución de los pesos intermedios estudio 1")
plot(density(pesosintermedios1), main="Densidad pesos intermedios estudio 1",
     xlab="Pesos intermedios", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(pesosintermedios1, main="Pesos intermedios estudio 1", outline=TRUE)

Peso medición

descr(Datos1[,6])
## Descriptive Statistics  
## Datos1$PesoFinal  
## N: 56  
## 
##                     PesoFinal
## ----------------- -----------
##              Mean      301.51
##           Std.Dev       35.85
##               Min      253.20
##                Q1      275.80
##            Median      294.75
##                Q3      314.20
##               Max      432.80
##               MAD       28.10
##               IQR       37.25
##                CV        0.12
##          Skewness        1.36
##       SE.Skewness        0.32
##          Kurtosis        2.12
##           N.Valid       56.00
##         Pct.Valid      100.00
pesosmedicion1<-Datos1[[6]]
par(mfrow = c(1, 2),cex=0.7)
hist(pesosmedicion1, col="pink", xlab="Pesos medición",
     ylab="Frecuencia", main="Distribución de los pesos medición estudio 1")
plot(density(pesosmedicion1), main="Densidad pesos medición estudio 1",
     xlab="Pesos medición", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(pesosmedicion1, main="Pesos medición estudio 1", outline=TRUE)

Tb. N

descr(Datos1[,8])
## Descriptive Statistics  
## Datos1$Tb.N  
## N: 56  
## 
##                       Tb.N
## ----------------- --------
##              Mean     3.14
##           Std.Dev     0.57
##               Min     1.86
##                Q1     2.76
##            Median     3.12
##                Q3     3.50
##               Max     4.62
##               MAD     0.54
##               IQR     0.73
##                CV     0.18
##          Skewness     0.08
##       SE.Skewness     0.32
##          Kurtosis    -0.27
##           N.Valid    56.00
##         Pct.Valid   100.00
TbN1<-na.omit(Datos1[[8]])
par(mfrow = c(1, 2),cex=0.7)
hist(TbN1, col="pink", xlab="Número trabecular",
     ylab="Frecuencia", main="Distribución del número trabecular estudio 1")
plot(density(TbN1), main="Densidad número trabecular estudio 1",
     xlab="Número trabecular", ylab="Densidad") 

par(mfrow = c(1, 1))
boxplot(TbN1, main="Número trabecular estudio 1", outline=TRUE)

BMC

descr(Datos1[,9])
## Descriptive Statistics  
## Datos1$BMC  
## N: 56  
## 
##                        BMC
## ----------------- --------
##              Mean     0.27
##           Std.Dev     0.03
##               Min     0.23
##                Q1     0.25
##            Median     0.27
##                Q3     0.29
##               Max     0.37
##               MAD     0.03
##               IQR     0.04
##                CV     0.12
##          Skewness     1.11
##       SE.Skewness     0.32
##          Kurtosis     0.67
##           N.Valid    56.00
##         Pct.Valid   100.00
BMC1<-na.omit(Datos1[[9]])
par(mfrow = c(1, 2),cex=0.7)
hist(BMC1, col="pink", xlab="BMC",
     ylab="Frecuencia", main="Distribución de BMC estudio 1")
plot(density(BMC1), main="Densidad BMC estudio 1",
     xlab="BMC", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(BMC1, main="BMC estudio 1", outline=TRUE)

BMD

descr(Datos1[,10])
## Descriptive Statistics  
## Datos1$BMD  
## N: 56  
## 
##                        BMD
## ----------------- --------
##              Mean   143.85
##           Std.Dev     7.46
##               Min   131.95
##                Q1   138.50
##            Median   142.71
##                Q3   147.78
##               Max   164.36
##               MAD     7.04
##               IQR     9.22
##                CV     0.05
##          Skewness     0.75
##       SE.Skewness     0.32
##          Kurtosis     0.34
##           N.Valid    56.00
##         Pct.Valid   100.00
BMD1<-na.omit(Datos1[[10]])
par(mfrow = c(1, 2),cex=0.7)
hist(BMD1, col="pink", xlab="BMD",
     ylab="Frecuencia", main="Distribución de BMD estudio 1")
plot(density(BMD1), main="BMD estudio 1",
     xlab="BMD", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(BMD1, main="BMD estudio 1", outline=TRUE)

Longitud

descr(Datos1[,11])
## Descriptive Statistics  
## Datos1$Longitud  
## N: 56  
## 
##                     Longitud
## ----------------- ----------
##              Mean      38.92
##           Std.Dev       1.10
##               Min      37.12
##                Q1      38.17
##            Median      38.74
##                Q3      39.51
##               Max      42.09
##               MAD       0.93
##               IQR       1.31
##                CV       0.03
##          Skewness       0.69
##       SE.Skewness       0.32
##          Kurtosis       0.00
##           N.Valid      56.00
##         Pct.Valid     100.00
longitud1<-na.omit(Datos1[[11]])
par(mfrow = c(1, 2),cex=0.7)
hist(longitud1, col="pink", xlab="Longitud",
     ylab="Frecuencia", main="Distribución de longitudes estudio 1")
plot(density(longitud1), main="Densidad de la longitud estudio 1",
     xlab="Longitud", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(longitud1, main="Longitudes estudio 1", outline=TRUE)

GPSurface

descr(Datos1[,12])
## Descriptive Statistics  
## Datos1$GPSurface  
## N: 56  
## 
##                     GPSurface
## ----------------- -----------
##              Mean       40.57
##           Std.Dev        3.07
##               Min       34.36
##                Q1       38.78
##            Median       40.34
##                Q3       42.62
##               Max       49.30
##               MAD        2.79
##               IQR        3.81
##                CV        0.08
##          Skewness        0.56
##       SE.Skewness        0.32
##          Kurtosis        0.35
##           N.Valid       56.00
##         Pct.Valid      100.00
GPSurface1<-na.omit(Datos1[[12]])
par(mfrow = c(1, 2),cex=0.7)
hist(GPSurface1, col="pink", xlab="GPSurface",
     ylab="Frecuencia", main="Distribución de GPSurface estudio 1")
plot(density(GPSurface1), main="Densidad GPSurface estudio 1",
     xlab="GPSurface", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(GPSurface1, main="GPSurface estudio 1", outline=TRUE)

GPVolume

descr(Datos1[,13])
## Descriptive Statistics  
## Datos1$GPVolume  
## N: 56  
## 
##                     GPVolume
## ----------------- ----------
##              Mean       5.60
##           Std.Dev       1.41
##               Min       3.35
##                Q1       4.66
##            Median       5.61
##                Q3       6.06
##               Max      10.32
##               MAD       1.17
##               IQR       1.38
##                CV       0.25
##          Skewness       1.06
##       SE.Skewness       0.32
##          Kurtosis       1.23
##           N.Valid      56.00
##         Pct.Valid     100.00
GPVolume1<-na.omit(Datos1[[13]])
par(mfrow = c(1, 2),cex=0.7)
hist(GPVolume1, col="pink", xlab="GPVolume",
     ylab="Frecuencia", main="Distribución de GPVolume estudio 1")
plot(density(GPVolume1), main="Densidad GPVolume estudio 1",
     xlab="GPVolume", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(GPVolume1, main="GPVolume estudio 1", outline=TRUE)

Tb.Sp

descr(Datos1[,7])
## Descriptive Statistics  
## Datos1$Tb.Sp  
## N: 56  
## 
##                      Tb.Sp
## ----------------- --------
##              Mean     0.33
##           Std.Dev     0.07
##               Min     0.21
##                Q1     0.28
##            Median     0.32
##                Q3     0.35
##               Max     0.55
##               MAD     0.05
##               IQR     0.07
##                CV     0.21
##          Skewness     0.93
##       SE.Skewness     0.32
##          Kurtosis     0.87
##           N.Valid    56.00
##         Pct.Valid   100.00
TbSp1<-na.omit(Datos1[[7]])
par(mfrow = c(1, 2),cex=0.7)
hist(TbSp1, col="pink", xlab="TbSp",
     ylab="Frecuencia", main="Distribución del espacio trabecular estudio 1")
plot(density(TbSp1), main="Densidad TbSp estudio 1",
     xlab="TbSp", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(TbSp1, main="TbSp estudio 1", outline=TRUE)

Outliers

boxplot(Datos1[,6:15],main="Outliers",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

Estandarizamos.

sca<-scale(Datos1[,6:15])

boxplot(sca,main="Outliers",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

ESTUDIO 2

Peso inicial

descr(Datos2[,14])
## Descriptive Statistics  
## Datos2$PesoInicial  
## N: 42  
## 
##                     PesoInicial
## ----------------- -------------
##              Mean        115.55
##           Std.Dev          8.18
##               Min         99.40
##                Q1        109.40
##            Median        116.55
##                Q3        119.90
##               Max        131.60
##               MAD          9.12
##               IQR         10.22
##                CV          0.07
##          Skewness          0.03
##       SE.Skewness          0.37
##          Kurtosis         -0.82
##           N.Valid         42.00
##         Pct.Valid        100.00
pesosiniciales2<-Datos2[[14]]
par(mfrow = c(1, 2),cex=0.7)
hist(pesosiniciales2, col="pink", xlab="Pesos iniciales", 
     ylab="Frecuencia", main="Distribución de los pesos iniciales estudio 2")
plot(density(pesosiniciales2), main="Densidad pesos iniciales estudio 2",
     xlab="Pesos iniciales", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(pesosiniciales2, main="Pesos iniciales estudio 2", outline=TRUE)

Peso intermedio

descr(Datos2[,15])
## Descriptive Statistics  
## Datos2$PesoIntermedio  
## N: 42  
## 
##                     PesoIntermedio
## ----------------- ----------------
##              Mean           201.07
##           Std.Dev            92.12
##               Min           123.20
##                Q1           144.90
##            Median           165.60
##                Q3           185.50
##               Max           397.80
##               MAD            29.73
##               IQR            39.35
##                CV             0.46
##          Skewness             1.38
##       SE.Skewness             0.37
##          Kurtosis             0.15
##           N.Valid            42.00
##         Pct.Valid           100.00
pesosintermedios2<-Datos2[[15]]
par(mfrow = c(1, 2),cex=0.7)
hist(pesosintermedios2, col="pink", xlab="Pesos intermedios", 
     ylab="Frecuencia", main="Distribución de los pesos intermedios estudio 2")
plot(density(pesosintermedios2), main="Densidad pesos intermedios estudio 2",
     xlab="Pesos intermedios", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(pesosintermedios2, main="Pesos intermedios estudio 2", outline=TRUE)

Peso medición

descr(Datos2[,6])
## Descriptive Statistics  
## Datos2$PesoFinal  
## N: 42  
## 
##                     PesoFinal
## ----------------- -----------
##              Mean      358.72
##           Std.Dev       54.66
##               Min      272.00
##                Q1      327.20
##            Median      345.30
##                Q3      370.60
##               Max      480.60
##               MAD       36.47
##               IQR       42.65
##                CV        0.15
##          Skewness        0.89
##       SE.Skewness        0.37
##          Kurtosis       -0.24
##           N.Valid       42.00
##         Pct.Valid      100.00
pesosmedicion2<-Datos2[[6]]
par(mfrow = c(1, 2),cex=0.7)
hist(pesosmedicion2, col="pink", xlab="Pesos medición",
     ylab="Frecuencia", main="Distribución de los pesos medición estudio 2")
plot(density(pesosmedicion2), main="Densidad pesos medición estudio 2",
     xlab="Pesos medición", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(pesosmedicion2, main="Pesos medición estudio 2", outline=TRUE)

Tb. N

descr(Datos2[,8])
## Descriptive Statistics  
## Datos2$Tb.N  
## N: 42  
## 
##                      Tb.N
## ----------------- -------
##              Mean    2.97
##           Std.Dev    0.61
##               Min    1.67
##                Q1    2.58
##            Median    3.02
##                Q3    3.34
##               Max    4.39
##               MAD    0.56
##               IQR    0.71
##                CV    0.20
##          Skewness   -0.09
##       SE.Skewness    0.37
##          Kurtosis   -0.28
##           N.Valid   40.00
##         Pct.Valid   95.24
TbN2<-na.omit(Datos2[[8]])
par(mfrow = c(1, 2),cex=0.7)
hist(TbN2, col="pink", xlab="Número trabecular",
     ylab="Frecuencia", main="Distribución del número trabecular estudio 2")
plot(density(TbN2), main="Densidad número trabecular estudio 2",
     xlab="Número trabecular", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(TbN2, main="Número trabecular estudio 2", outline=TRUE)

BMC

descr(Datos2[,9])
## Descriptive Statistics  
## Datos2$BMC  
## N: 42  
## 
##                       BMC
## ----------------- -------
##              Mean    0.35
##           Std.Dev    0.05
##               Min    0.26
##                Q1    0.32
##            Median    0.33
##                Q3    0.35
##               Max    0.50
##               MAD    0.03
##               IQR    0.03
##                CV    0.15
##          Skewness    1.19
##       SE.Skewness    0.38
##          Kurtosis    0.72
##           N.Valid   38.00
##         Pct.Valid   90.48
BMC2<-na.omit(Datos2[[9]])
par(mfrow = c(1, 2),cex=0.7)
hist(BMC2, col="pink", xlab="BMC",
     ylab="Frecuencia", main="Distribución de BMC estudio 2")
plot(density(BMC2), main="Densidad BMC estudio 2",
     xlab="BMC", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(BMC2, main="BMC estudio 2", outline=TRUE)

BMD

descr(Datos2[,10])
## Descriptive Statistics  
## Datos2$BMD  
## N: 42  
## 
##                        BMD
## ----------------- --------
##              Mean   159.34
##           Std.Dev    12.07
##               Min   137.87
##                Q1   151.25
##            Median   157.38
##                Q3   166.30
##               Max   186.94
##               MAD     9.45
##               IQR    14.53
##                CV     0.08
##          Skewness     0.66
##       SE.Skewness     0.38
##          Kurtosis    -0.36
##           N.Valid    38.00
##         Pct.Valid    90.48
BMD2<-na.omit(Datos2[[10]])
par(mfrow = c(1, 2),cex=0.7)
hist(BMD2, col="pink", xlab="BMD",
     ylab="Frecuencia", main="Distribución de BMD estudio 2")
plot(density(BMD2), main="BMD estudio 2",
     xlab="BMD", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(BMD2, main="BMD estudio 2", outline=TRUE)

Longitud

descr(Datos2[,11])
## Descriptive Statistics  
## Datos2$Longitud  
## N: 42  
## 
##                     Longitud
## ----------------- ----------
##              Mean      42.04
##           Std.Dev       1.47
##               Min      39.29
##                Q1      41.32
##            Median      41.74
##                Q3      42.73
##               Max      45.52
##               MAD       0.73
##               IQR       1.27
##                CV       0.04
##          Skewness       0.65
##       SE.Skewness       0.38
##          Kurtosis       0.04
##           N.Valid      38.00
##         Pct.Valid      90.48
longitud2<-na.omit(Datos2[[11]])
par(mfrow = c(1, 2),cex=0.7)
hist(longitud2, col="pink", xlab="Longitud",
     ylab="Frecuencia", main="Distribución de longitudes estudio 2")
plot(density(longitud2), main="Densidad de la longitud estudio 2",
     xlab="Longitud", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(longitud2, main="Longitudes estudio 2", outline=TRUE)

GPSurface

descr(Datos2[,12])
## Descriptive Statistics  
## Datos2$GPSurface  
## N: 42  
## 
##                     GPSurface
## ----------------- -----------
##              Mean       44.12
##           Std.Dev        4.44
##               Min       35.71
##                Q1       41.27
##            Median       44.03
##                Q3       45.63
##               Max       54.68
##               MAD        3.62
##               IQR        4.28
##                CV        0.10
##          Skewness        0.46
##       SE.Skewness        0.37
##          Kurtosis       -0.15
##           N.Valid       40.00
##         Pct.Valid       95.24
GPSurface2<-na.omit(Datos2[[12]])
par(mfrow = c(1, 2),cex=0.7)
hist(GPSurface2, col="pink", xlab="GPSurface",
     ylab="Frecuencia", main="Distribución de GPSurface estudio 2")
plot(density(GPSurface2), main="Densidad GPSurface estudio 2",
     xlab="GPSurface", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(GPSurface2, main="GPSurface estudio 2", outline=TRUE)

GPVolume

descr(Datos2[,13])
## Descriptive Statistics  
## Datos2$GPVolume  
## N: 42  
## 
##                     GPVolume
## ----------------- ----------
##              Mean       5.57
##           Std.Dev       0.97
##               Min       3.59
##                Q1       4.79
##            Median       5.52
##                Q3       6.15
##               Max       7.85
##               MAD       0.99
##               IQR       1.33
##                CV       0.17
##          Skewness       0.26
##       SE.Skewness       0.37
##          Kurtosis      -0.41
##           N.Valid      40.00
##         Pct.Valid      95.24
GPVolume2<-na.omit(Datos2[[13]])
par(mfrow = c(1, 2),cex=0.7)
hist(GPVolume2, col="pink", xlab="GPVolume",
     ylab="Frecuencia", main="Distribución de GPVolume estudio 2")
plot(density(GPVolume2), main="Densidad GPVolume estudio 2",
     xlab="GPVolume", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(GPVolume2, main="GPVolume estudio 2", outline=TRUE)

Tb.Sp

descr(Datos2[,7])
## Descriptive Statistics  
## Datos2$Tb.Sp  
## N: 42  
## 
##                     Tb.Sp
## ----------------- -------
##              Mean    0.35
##           Std.Dev    0.09
##               Min    0.21
##                Q1    0.29
##            Median    0.33
##                Q3    0.38
##               Max    0.60
##               MAD    0.06
##               IQR    0.08
##                CV    0.25
##          Skewness    1.09
##       SE.Skewness    0.37
##          Kurtosis    0.80
##           N.Valid   40.00
##         Pct.Valid   95.24
TbSp2<-na.omit(Datos2[[7]])
par(mfrow = c(1, 2),cex=0.7)
hist(TbSp2, col="pink", xlab="TbSp",
     ylab="Frecuencia", main="Distribución del espacio trabecular estudio 2")
plot(density(TbSp2), main="Densidad TbSp estudio 2",
     xlab="TbSp", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(TbSp2, main="TbSp estudio 2", outline=TRUE)

Outliers

boxplot(Datos2[,6:15],main="Outliers",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

Estandarizamos.

sca<-scale(Datos2[,6:15])

boxplot(sca,main="Outliers",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

ESTUDIO 3

Peso inicial

descr(Datos3[,14])
## Descriptive Statistics  
## Datos3$PesoInicial  
## N: 36  
## 
##                     PesoInicial
## ----------------- -------------
##              Mean         71.76
##           Std.Dev          4.56
##               Min         63.00
##                Q1         69.35
##            Median         71.65
##                Q3         73.30
##               Max         86.80
##               MAD          2.74
##               IQR          3.83
##                CV          0.06
##          Skewness          0.95
##       SE.Skewness          0.39
##          Kurtosis          2.16
##           N.Valid         36.00
##         Pct.Valid        100.00
pesosiniciales3<-Datos3[[14]]
par(mfrow = c(1, 2),cex=0.7)
hist(pesosiniciales3, col="pink", xlab="Pesos iniciales", 
     ylab="Frecuencia", main="Distribución de los pesos iniciales estudio 3")
plot(density(pesosiniciales3), main="Densidad pesos iniciales estudio 3",
     xlab="Pesos iniciales", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(pesosiniciales3, main="Pesos iniciales estudio 3", outline=TRUE)

Peso intermedio

descr(Datos3[,15])
## Descriptive Statistics  
## Datos3$PesoIntermedio  
## N: 36  
## 
##                     PesoIntermedio
## ----------------- ----------------
##              Mean           127.51
##           Std.Dev            51.49
##               Min            77.00
##                Q1            95.05
##            Median           104.25
##                Q3           166.65
##               Max           234.00
##               MAD            17.27
##               IQR            54.10
##                CV             0.40
##          Skewness             1.04
##       SE.Skewness             0.39
##          Kurtosis            -0.73
##           N.Valid            36.00
##         Pct.Valid           100.00
pesosintermedios3<-Datos3[[15]]
par(mfrow = c(1, 2),cex=0.7)
hist(pesosintermedios3, col="pink", xlab="Pesos intermedios", 
     ylab="Frecuencia", main="Distribución de los pesos intermedios estudio 3")
plot(density(pesosintermedios3), main="Densidad pesos intermedios estudio 3",
     xlab="Pesos intermedios", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(pesosintermedios3, main="Pesos intermedios estudio 3", outline=TRUE)

Peso medición

descr(Datos3[,6])
## Descriptive Statistics  
## Datos3$PesoFinal  
## N: 36  
## 
##                     PesoFinal
## ----------------- -----------
##              Mean      238.28
##           Std.Dev       47.10
##               Min      163.00
##                Q1      208.00
##            Median      223.00
##                Q3      275.65
##               Max      330.00
##               MAD       26.02
##               IQR       61.32
##                CV        0.20
##          Skewness        0.68
##       SE.Skewness        0.39
##          Kurtosis       -0.84
##           N.Valid       36.00
##         Pct.Valid      100.00
pesosmedicion3<-Datos3[[6]]
par(mfrow = c(1, 2),cex=0.7)
hist(pesosmedicion3, col="pink", xlab="Pesos medición",
     ylab="Frecuencia", main="Distribución de los pesos medición estudio 3")
plot(density(pesosmedicion3), main="Densidad pesos medición estudio 3",
     xlab="Pesos medición", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(pesosmedicion3, main="Pesos medición estudio 3", outline=TRUE)

Tb. N

descr(Datos3[,8])
## Descriptive Statistics  
## Datos3$Tb.N  
## N: 36  
## 
##                       Tb.N
## ----------------- --------
##              Mean     3.29
##           Std.Dev     0.61
##               Min     1.79
##                Q1     2.78
##            Median     3.48
##                Q3     3.75
##               Max     4.24
##               MAD     0.69
##               IQR     0.95
##                CV     0.18
##          Skewness    -0.45
##       SE.Skewness     0.39
##          Kurtosis    -0.73
##           N.Valid    36.00
##         Pct.Valid   100.00
TbN3<-na.omit(Datos3[[8]])
par(mfrow = c(1, 2),cex=0.7)
hist(TbN3, col="pink", xlab="Número trabecular",
     ylab="Frecuencia", main="Distribución del número trabecular estudio 3")
plot(density(TbN3), main="Densidad número trabecular estudio 3",
     xlab="Número trabecular", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(TbN3, main="Número trabecular estudio 3", outline=TRUE)

BMC

descr(Datos3[,9])
## Descriptive Statistics  
## Datos3$BMC  
## N: 36  
## 
##                        BMC
## ----------------- --------
##              Mean     0.20
##           Std.Dev     0.04
##               Min     0.15
##                Q1     0.17
##            Median     0.19
##                Q3     0.24
##               Max     0.30
##               MAD     0.03
##               IQR     0.06
##                CV     0.21
##          Skewness     0.76
##       SE.Skewness     0.39
##          Kurtosis    -0.73
##           N.Valid    36.00
##         Pct.Valid   100.00
BMC3<-na.omit(Datos3[[9]])
par(mfrow = c(1, 2),cex=0.7)
hist(BMC3, col="pink", xlab="BMC",
     ylab="Frecuencia", main="Distribución de BMC estudio 3")
plot(density(BMC3), main="Densidad BMC estudio 3",
     xlab="BMC", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(BMC3, main="BMC estudio 3", outline=TRUE)

BMD

descr(Datos3[,10])
## Descriptive Statistics  
## Datos3$BMD  
## N: 36  
## 
##                        BMD
## ----------------- --------
##              Mean   123.57
##           Std.Dev    11.05
##               Min   105.62
##                Q1   113.55
##            Median   122.86
##                Q3   131.92
##               Max   146.00
##               MAD    13.56
##               IQR    18.03
##                CV     0.09
##          Skewness     0.42
##       SE.Skewness     0.39
##          Kurtosis    -0.93
##           N.Valid    36.00
##         Pct.Valid   100.00
BMD3<-na.omit(Datos3[[10]])
par(mfrow = c(1, 2),cex=0.7)
hist(BMD3, col="pink", xlab="BMD",
     ylab="Frecuencia", main="Distribución de BMD estudio 3")
plot(density(BMD3), main="BMD estudio 3",
     xlab="BMD", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(BMD3, main="BMD estudio 3", outline=TRUE)

Longitud

descr(Datos3[,11])
## Descriptive Statistics  
## Datos3$Longitud  
## N: 36  
## 
##                     Longitud
## ----------------- ----------
##              Mean      34.84
##           Std.Dev       1.59
##               Min      32.45
##                Q1      33.61
##            Median      34.30
##                Q3      36.16
##               Max      37.89
##               MAD       1.19
##               IQR       2.46
##                CV       0.05
##          Skewness       0.55
##       SE.Skewness       0.39
##          Kurtosis      -1.02
##           N.Valid      36.00
##         Pct.Valid     100.00
longitud3<-na.omit(Datos3[[11]])
par(mfrow = c(1, 2),cex=0.7)
hist(longitud3, col="pink", xlab="Longitud",
     ylab="Frecuencia", main="Distribución de longitudes estudio 3")
plot(density(longitud3), main="Densidad de la longitud estudio 3",
     xlab="Longitud", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(longitud3, main="Longitudes estudio 3", outline=TRUE)

GPSurface

descr(Datos3[,12])
## Descriptive Statistics  
## Datos3$GPSurface  
## N: 36  
## 
##                     GPSurface
## ----------------- -----------
##              Mean       35.90
##           Std.Dev        3.71
##               Min       28.96
##                Q1       33.73
##            Median       34.47
##                Q3       38.31
##               Max       44.20
##               MAD        1.93
##               IQR        4.43
##                CV        0.10
##          Skewness        0.69
##       SE.Skewness        0.39
##          Kurtosis       -0.36
##           N.Valid       36.00
##         Pct.Valid      100.00
GPSurface3<-na.omit(Datos3[[12]])
par(mfrow = c(1, 2),cex=0.7)
hist(GPSurface3, col="pink", xlab="GPSurface",
     ylab="Frecuencia", main="Distribución de GPSurface estudio 3")
plot(density(GPSurface3), main="Densidad GPSurface estudio 3",
     xlab="GPSurface", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(GPSurface3, main="GPSurface estudio 3", outline=TRUE)

GPVolume

descr(Datos3[,13])
## Descriptive Statistics  
## Datos3$GPVolume  
## N: 36  
## 
##                     GPVolume
## ----------------- ----------
##              Mean       5.65
##           Std.Dev       0.86
##               Min       4.04
##                Q1       4.99
##            Median       5.57
##                Q3       6.32
##               Max       7.58
##               MAD       0.93
##               IQR       1.28
##                CV       0.15
##          Skewness       0.13
##       SE.Skewness       0.39
##          Kurtosis      -0.69
##           N.Valid      36.00
##         Pct.Valid     100.00
GPVolume3<-na.omit(Datos3[[13]])
par(mfrow = c(1, 2),cex=0.7)
hist(GPVolume3, col="pink", xlab="GPVolume",
     ylab="Frecuencia", main="Distribución de GPVolume estudio 3")
plot(density(GPVolume3), main="Densidad GPVolume estudio 3",
     xlab="GPVolume", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(GPVolume3, main="GPVolume estudio 3", outline=TRUE)

Tb.Sp

descr(Datos3[,7])
## Descriptive Statistics  
## Datos3$Tb.Sp  
## N: 36  
## 
##                      Tb.Sp
## ----------------- --------
##              Mean     0.31
##           Std.Dev     0.07
##               Min     0.22
##                Q1     0.26
##            Median     0.28
##                Q3     0.36
##               Max     0.57
##               MAD     0.05
##               IQR     0.10
##                CV     0.24
##          Skewness     1.25
##       SE.Skewness     0.39
##          Kurtosis     1.57
##           N.Valid    36.00
##         Pct.Valid   100.00
TbSp3<-na.omit(Datos3[[7]])
par(mfrow = c(1, 2),cex=0.7)
hist(TbSp3, col="pink", xlab="TbSp",
     ylab="Frecuencia", main="Distribución del espacio trabecular estudio 3")
plot(density(TbSp3), main="Densidad TbSp estudio 3",
     xlab="TbSp", ylab="Densidad")

par(mfrow = c(1, 1))
boxplot(TbSp3, main="TbSp estudio 3", outline=TRUE)

Outliers

boxplot(Datos3[,6:15],main="Outliers",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

Estandarizamos.

sca<-scale(Datos3[,6:15])

boxplot(sca,main="Outliers",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

Vamos a superponer los histogramas para ver cómo cambian las variables en función de la dosis.

Estudio 1

par(mfrow = c(1, 2),cex=0.7)
p1_1<-ggplot(data=Datos1,aes(x=PesoInicial,fill=factor(Dosis))) +
    geom_histogram(position = "identity", alpha = 0.5, bins=20)
p2_1<-ggplot(data=Datos1,aes(x=PesoIntermedio,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p3_1<-ggplot(data=Datos1,aes(x=PesoFinal,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p4_1<-ggplot(data=Datos1,aes(x=Tb.N,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p5_1<-ggplot(data=Datos1,aes(x=Tb.Sp,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p6_1<-ggplot(data=Datos1,aes(x=BMC,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p7_1<-ggplot(data=Datos1,aes(x=BMD,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p8_1<-ggplot(data=Datos1,aes(x=Longitud,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p9_1<-ggplot(data=Datos1,aes(x=GPSurface,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p10_1<-ggplot(data=Datos1,aes(x=GPVolume,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)

ggarrange(p1_1,p2_1, p3_1,p4_1,p5_1,p6_1,p7_1,p8_1,p9_1,p10_1, nrow = 4, ncol=3, common.legend = TRUE)

Estudio 2

Datos2.2<-na.omit(Datos2)
par(mfrow = c(1, 2),cex=0.7)
p1_2<-ggplot(data=Datos2.2,aes(x=PesoInicial,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p2_2<-ggplot(data=Datos2.2,aes(x=PesoIntermedio,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p3_2<-ggplot(data=Datos2.2,aes(x=PesoFinal,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p4_2<-ggplot(data=Datos2.2,aes(x=Tb.N,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p5_2<-ggplot(data=Datos2.2,aes(x=Tb.Sp,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p6_2<-ggplot(data=Datos2.2,aes(x=BMC,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p7_2<-ggplot(data=Datos2.2,aes(x=BMD,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p8_2<-ggplot(data=Datos2.2,aes(x=Longitud,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p9_2<-ggplot(data=Datos2.2,aes(x=GPSurface,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p10_2<-ggplot(data=Datos2.2,aes(x=GPVolume,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)

ggarrange(p1_2,p2_2, p3_2,p4_2,p5_2,p6_2,p7_2,p8_2,p9_2,p10_2, nrow = 4, ncol=3, common.legend = TRUE)

Estudio 3

par(mfrow = c(1, 2),cex=0.7)
p1_3<-ggplot(data=Datos3,aes(x=PesoInicial,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p2_3<-ggplot(data=Datos3,aes(x=PesoIntermedio,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p3_3<-ggplot(data=Datos3,aes(x=PesoFinal,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p4_3<-ggplot(data=Datos3,aes(x=Tb.N,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p5_3<-ggplot(data=Datos3,aes(x=Tb.Sp,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p6_3<-ggplot(data=Datos3,aes(x=BMC,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p7_3<-ggplot(data=Datos3,aes(x=BMD,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p8_3<-ggplot(data=Datos3,aes(x=Longitud,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p9_3<-ggplot(data=Datos3,aes(x=GPSurface,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)
p10_3<-ggplot(data=Datos3,aes(x=GPVolume,fill=factor(Dosis))) +
  geom_histogram(position = "identity", alpha = 0.5, bins=20)

ggarrange(p1_3,p2_3, p3_3,p4_3,p5_3,p6_3,p7_3,p8_3,p9_3,p10_3, nrow = 4, ncol=3, common.legend = TRUE)

En general parece que el peso intermedio, el BMC y la longitud separan mejor al grupo de control. Ahora vamos a realizar un análisis exploratorio visual bivariante.

Análisis exploratorio bivariante

Estudio 1

pairs(Datos1[, c("PesoInicial", "PesoIntermedio", "PesoFinal", "Tb.N", "Tb.Sp", "BMC", "BMD", "Longitud", "GPSurface", "GPVolume")],
      col = rainbow(length(levels(factor(Datos1$Dosis))))[factor(Datos1$Dosis)],
      pch = 19)
legend("topright", 
       legend = levels(factor(Datos1$Dosis)), 
       fill = rainbow(length(levels(factor(Datos1$Dosis)))))

Correlaciones

par(mar = c(3, 1, 4, 1))
corrplot(cor(Datos1[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Estudio 1", side = 3, line = 3, adj = 0.3, cex = 1.1)

Correlaciones por dosis

dato1dosis0<-subset(Datos1,Dosis==0)
dato1dosis1<-subset(Datos1,Dosis==1)
dato1dosis8<-subset(Datos1,Dosis==8)
dato1dosis10<-subset(Datos1,Dosis==10)
dato1dosis80<-subset(Datos1,Dosis==80)
dato1dosis100<-subset(Datos1,Dosis==100)

par(mar = c(3, 1, 4, 1))
corrplot(cor(dato1dosis0[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Dosis 0", side = 3, line = 3, adj = 0.3, cex = 1.1)

corrplot(cor(dato1dosis1[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Dosis 1", side = 3, line = 3, adj = 0.3, cex = 1.1)

corrplot(cor(dato1dosis8[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Dosis 8", side = 3, line = 3, adj = 0.3, cex = 1.1)

corrplot(cor(dato1dosis10[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Dosis 10", side = 3, line = 3, adj = 0.3, cex = 1.1)

corrplot(cor(dato1dosis80[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Dosis 80", side = 3, line = 3, adj = 0.3, cex = 1.1)

corrplot(cor(dato1dosis100[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Dosis 100", side = 3, line = 3, adj = 0.3, cex = 1.1)

Estudio 2

Omitimos los valores perdidos del Estudio 2 para hacer la gráfica. No sería un error eliminarlos ya que son menos del 5%.

Datos2.2<-na.omit(Datos2)
pairs(Datos2.2[, c("PesoInicial", "PesoIntermedio", "PesoFinal", "Tb.N", "Tb.Sp", "BMC", "BMD", "Longitud", "GPSurface", "GPVolume")],
      col = rainbow(length(levels(factor(Datos2.2$Dosis))))[factor(Datos2.2$Dosis)],pch = 19)
legend("topright", 
       legend = levels(factor(Datos2.2$Dosis)), 
       fill = rainbow(length(levels(factor(Datos2.2$Dosis)))))

Correlaciones

par(mar = c(3, 1, 4, 1))
corrplot(cor(Datos2.2[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Estudio 2", side = 3, line = 3, adj = 0.3, cex = 1.1)

Correlaciones por dosis

dato2dosis0<-subset(Datos2.2,Dosis==0)
dato2dosis80<-subset(Datos2.2,Dosis==80)
dato2dosis400<-subset(Datos2.2,Dosis==400)

par(mar = c(3, 1, 4, 1))
corrplot(cor(dato2dosis0[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Dosis 0", side = 3, line = 3, adj = 0.3, cex = 1.1)

corrplot(cor(dato2dosis80[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Dosis 80", side = 3, line = 3, adj = 0.3, cex = 1.1)

corrplot(cor(dato2dosis400[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Dosis 400", side = 3, line = 3, adj = 0.3, cex = 1.1)

Estudio 3

pairs(Datos3[, c("PesoInicial", "PesoIntermedio", "PesoFinal", "Tb.N", "Tb.Sp", "BMC", "BMD", "Longitud", "GPSurface", "GPVolume")],
      col = rainbow(length(levels(factor(Datos3$Dosis))))[factor(Datos3$Dosis)],pch = 19)
legend("topright", 
       legend = levels(factor(Datos3$Dosis)), 
       fill = rainbow(length(levels(factor(Datos3$Dosis)))))

#### Correlaciones

par(mar = c(3, 1, 4, 1))
corrplot(cor(Datos3[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Estudio 3", side = 3, line = 2, adj = 0.3, cex = 1.1)

Correlaciones por dosis

dato3dosis0<-subset(Datos3,Dosis==0)
dato3dosis500<-subset(Datos3,Dosis==500)
dato3dosis1000<-subset(Datos3,Dosis==1000)

par(mar = c(3, 1, 4, 1))
corrplot(cor(dato3dosis0[6:15]), order = "hclust", tl.col='black', tl.cex=0.5)
mtext("Dosis 0", side = 3, line = 3, adj = 0.3, cex = 1.1)

corrplot(cor(dato3dosis500[6:15]), order = "hclust", tl.col='black', tl.cex=0.5)
mtext("Dosis 500", side = 3, line = 3, adj = 0.3, cex = 1.1)

corrplot(cor(dato3dosis1000[6:15]), order = "hclust", tl.col='black', tl.cex=0.5)
mtext("Dosis 1000", side = 3, line = 3, adj = 0.3, cex = 1.1)

Análisis exploratorio trivariante

Finalmente, nos queda hacer un análisis exploratorio visual trivariante. Comparamos peso intermedio, longitud, BMC, BMD, GPSurface y Tb.Sp que en los gráficos de histogramas parecían los que más información nos aportaban.

Estudio 1

Datos1Factor <- Datos1
Datos2.2Factor <- Datos2.2
Datos3Factor <- Datos3
Datos1Factor$Dosis <- factor(Datos1Factor$Dosis)
Datos2.2Factor$Dosis <- factor(Datos2.2Factor$Dosis)
Datos3Factor$Dosis <- factor(Datos3Factor$Dosis)

par(mfrow = c(1, 2),cex=0.7)
#####
scatterplot3d(Datos1$PesoInicial, Datos1$PesoIntermedio, Datos1$PesoFinal,
              color = rainbow(length(levels(Datos1Factor$Dosis)))[as.factor(Datos1$Dosis)],
              pch = 19, grid = TRUE, xlab = "Peso Inicial", ylab = "Peso Intermedio", zlab = "Peso Final",
              angle = 65, cex.axis = 0.6)

legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(factor(Datos1$Dosis)), fill = rainbow(length(levels(factor(Datos1$Dosis)))))
#####
scatterplot3d(Datos1$Longitud, Datos1$PesoIntermedio, Datos1$BMC,
              color = rainbow(length(levels(Datos1Factor$Dosis)))[as.factor(Datos1$Dosis)], pch = 19,
              grid = TRUE, xlab = "Longitud", ylab = "Peso Intermedio", zlab = "BMC",
              angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(Datos1Factor$Dosis), fill = rainbow(length(levels(Datos1Factor$Dosis))))

####
scatterplot3d(Datos1Factor$GPSurface, Datos1Factor$PesoIntermedio, Datos1Factor$BMC,
              color = rainbow(length(levels(Datos1Factor$Dosis)))[Datos1Factor$Dosis], pch = 19,
              grid = TRUE, xlab = "GPSurface", ylab = "Peso Intermedio", zlab = "BMC",
              angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(Datos1Factor$Dosis), fill = rainbow(length(levels(Datos1Factor$Dosis))))
####
scatterplot3d(Datos1Factor$GPSurface, Datos1Factor$PesoIntermedio, Datos1Factor$Longitud,
              color = rainbow(length(levels(Datos1Factor$Dosis)))[Datos1Factor$Dosis], pch = 19,
              grid = TRUE, xlab = "GPSurface", ylab = "Peso Intermedio", zlab = "Longitud",
              angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(Datos1Factor$Dosis), fill = rainbow(length(levels(Datos1Factor$Dosis))))

####
scatterplot3d(Datos1Factor$BMD, Datos1Factor$PesoIntermedio, Datos1Factor$Tb.Sp,
              color = rainbow(length(levels(Datos1Factor$Dosis)))[Datos1Factor$Dosis], pch = 19,
              grid = TRUE, xlab = "BMD", ylab = "Peso Intermedio", zlab = "Tb.Sp",
              angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(Datos1Factor$Dosis), fill = rainbow(length(levels(Datos1Factor$Dosis))))
#####
scatterplot3d(Datos1Factor$BMD, Datos1Factor$PesoIntermedio, Datos1Factor$Longitud,
              color = rainbow(length(levels(Datos1Factor$Dosis)))[Datos1Factor$Dosis], pch = 19,
              grid = TRUE, xlab = "BMD", ylab = "Peso Intermedio", zlab = "Longitud",
              angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(Datos1Factor$Dosis), fill = rainbow(length(levels(Datos1Factor$Dosis))))

#####
scatterplot3d(Datos1Factor$BMC, Datos1Factor$PesoIntermedio, Datos1Factor$Tb.Sp,
              color = rainbow(length(levels(Datos1Factor$Dosis)))[Datos1Factor$Dosis], pch = 19,
              grid = TRUE, xlab = "BMC", ylab = "Peso Intermedio", zlab = "Tb.Sp",
              angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(Datos1Factor$Dosis), fill = rainbow(length(levels(Datos1Factor$Dosis))))
#####
scatterplot3d(Datos1Factor$BMD, Datos1Factor$PesoIntermedio, Datos1Factor$GPSurface,
              color = rainbow(length(levels(Datos1Factor$Dosis)))[Datos1Factor$Dosis], pch = 19,
              grid = TRUE, xlab = "BMD", ylab = "Peso Intermedio", zlab = "GPSurface",
              angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(Datos1Factor$Dosis), fill = rainbow(length(levels(Datos1Factor$Dosis))))

Estudio 2

par(mfrow = c(1, 2),cex=0.7)

scatterplot3d(Datos2.2$PesoInicial, Datos2.2$PesoIntermedio, Datos2.2$PesoFinal,
              color = rainbow(length(levels(factor(Datos2.2$Dosis))))[as.factor(Datos2.2$Dosis)],
              pch = 19, grid = TRUE, xlab = "Peso Inicial", ylab = "Peso Intermedio", zlab = "Peso Final",
              angle = 65, cex.axis = 0.6)

legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(factor(Datos2.2$Dosis)), fill = rainbow(length(levels(factor(Datos2.2$Dosis)))))
#####
scatterplot3d(Datos2.2Factor$Longitud, Datos2.2Factor$PesoIntermedio, Datos2.2Factor$BMC,
              color = rainbow(length(levels(Datos2.2Factor$Dosis)))[Datos2.2Factor$Dosis], pch = 19,
              grid = TRUE, xlab = "Longitud", ylab = "Peso Intermedio", zlab = "BMC",
              angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(Datos2.2Factor$Dosis), fill = rainbow(length(levels(Datos2.2Factor$Dosis))))

#####
scatterplot3d(Datos2.2Factor$GPSurface, Datos2.2Factor$PesoIntermedio, Datos2.2Factor$BMC,
              color = rainbow(length(levels(Datos2.2Factor$Dosis)))[Datos2.2Factor$Dosis], pch = 19,
              grid = TRUE, xlab = "GPSurface", ylab = "Peso Intermedio", zlab = "BMC",
              angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(Datos2.2Factor$Dosis), fill = rainbow(length(levels(Datos2.2Factor$Dosis))))
######
scatterplot3d(Datos2.2Factor$GPSurface, Datos2.2Factor$PesoIntermedio, Datos2.2Factor$Longitud,
              color = rainbow(length(levels(Datos2.2Factor$Dosis)))[Datos2.2Factor$Dosis], pch = 19,
              grid = TRUE, xlab = "GPSurface", ylab = "Peso Intermedio", zlab = "Longitud",
              angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(Datos2.2Factor$Dosis), fill = rainbow(length(levels(Datos2.2Factor$Dosis))))

#######
scatterplot3d(Datos2.2Factor$BMD, Datos2.2Factor$PesoIntermedio, Datos2.2Factor$Tb.Sp,
              color = rainbow(length(levels(Datos2.2Factor$Dosis)))[Datos2.2Factor$Dosis], pch = 19,
              grid = TRUE, xlab = "BMD", ylab = "Peso Intermedio", zlab = "Tb.Sp",
              angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(Datos2.2Factor$Dosis), fill = rainbow(length(levels(Datos2.2Factor$Dosis))))
######
scatterplot3d(Datos2.2Factor$BMD, Datos2.2Factor$PesoIntermedio, Datos2.2Factor$Longitud,
              color = rainbow(length(levels(Datos2.2Factor$Dosis)))[Datos2.2Factor$Dosis], pch = 19,
              grid = TRUE, xlab = "BMD", ylab = "Peso Intermedio", zlab = "Longitud",
              angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(Datos2.2Factor$Dosis), fill = rainbow(length(levels(Datos2.2Factor$Dosis))))

#######
scatterplot3d(Datos2.2Factor$BMC, Datos2.2Factor$PesoIntermedio, Datos2.2Factor$Tb.Sp,
              color = rainbow(length(levels(Datos2.2Factor$Dosis)))[Datos2.2Factor$Dosis], pch = 19,
              grid = TRUE, xlab = "BMC", ylab = "Peso Intermedio", zlab = "Tb.Sp",
              angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(Datos2.2Factor$Dosis), fill = rainbow(length(levels(Datos2.2Factor$Dosis))))
######
scatterplot3d(Datos2.2Factor$BMD, Datos2.2Factor$PesoIntermedio, Datos2.2Factor$GPSurface,
              color = rainbow(length(levels(Datos2.2Factor$Dosis)))[Datos2.2Factor$Dosis], pch = 19,
              grid = TRUE, xlab = "BMD", ylab = "Peso Intermedio", zlab = "GPSurface",
              angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(Datos2.2Factor$Dosis), fill = rainbow(length(levels(Datos2.2Factor$Dosis))))

Estudio 3

par(mfrow = c(1, 2),cex=0.7)

scatterplot3d(Datos3$PesoInicial, Datos3$PesoIntermedio, Datos3$PesoFinal,
              color = rainbow(length(levels(factor(Datos3$Dosis))))[as.factor(Datos3$Dosis)],
              pch = 19, grid = TRUE, xlab = "Peso Inicial", ylab = "Peso Intermedio", zlab = "Peso Final",
              angle = 65, cex.axis = 0.6)

legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(factor(Datos3$Dosis)), fill = rainbow(length(levels(factor(Datos3$Dosis)))))
#####
scatterplot3d(Datos3Factor$Longitud, Datos3Factor$PesoIntermedio, Datos3Factor$BMC,
              color = rainbow(length(levels(Datos3Factor$Dosis)))[Datos3Factor$Dosis], pch = 19,
              grid = TRUE, xlab = "Longitud", ylab = "Peso Intermedio", zlab = "BMC",
              angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(Datos3Factor$Dosis), fill = rainbow(length(levels(Datos3Factor$Dosis))))

######
scatterplot3d(Datos3Factor$GPSurface, Datos3Factor$PesoIntermedio, Datos3Factor$BMC,
              color = rainbow(length(levels(Datos3Factor$Dosis)))[Datos3Factor$Dosis], pch = 19,
              grid = TRUE, xlab = "GPSurface", ylab = "Peso Intermedio", zlab = "BMC",
              angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(Datos3Factor$Dosis), fill = rainbow(length(levels(Datos3Factor$Dosis))))
#######
scatterplot3d(Datos3Factor$GPSurface, Datos3Factor$PesoIntermedio, Datos3Factor$Longitud,
              color = rainbow(length(levels(Datos3Factor$Dosis)))[Datos3Factor$Dosis], pch = 19,
              grid = TRUE, xlab = "GPSurface", ylab = "Peso Intermedio", zlab = "Longitud",
              angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(Datos3Factor$Dosis), fill = rainbow(length(levels(Datos3Factor$Dosis))))

#########
scatterplot3d(Datos3Factor$BMD, Datos3Factor$PesoIntermedio, Datos3Factor$Tb.Sp,
              color = rainbow(length(levels(Datos3Factor$Dosis)))[Datos3Factor$Dosis], pch = 19,
              grid = TRUE, xlab = "GPSurface", ylab = "Peso Intermedio", zlab = "Tb.Sp",
              angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(Datos3Factor$Dosis), fill = rainbow(length(levels(Datos3Factor$Dosis))))
#########
scatterplot3d(Datos3Factor$BMD, Datos3Factor$PesoIntermedio, Datos3Factor$Longitud,
              color = rainbow(length(levels(Datos3Factor$Dosis)))[Datos3Factor$Dosis], pch = 19,
              grid = TRUE, xlab = "BMD", ylab = "Peso Intermedio", zlab = "Longitud",
              angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(Datos3Factor$Dosis), fill = rainbow(length(levels(Datos3Factor$Dosis))))

#######
scatterplot3d(Datos3Factor$BMD, Datos3Factor$PesoIntermedio, Datos3Factor$Tb.Sp,
              color = rainbow(length(levels(Datos3Factor$Dosis)))[Datos3Factor$Dosis], pch = 19,
              grid = TRUE, xlab = "GPSurface", ylab = "Peso Intermedio", zlab = "Tb.Sp",
              angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(Datos3Factor$Dosis), fill = rainbow(length(levels(Datos3Factor$Dosis))))
#######
scatterplot3d(Datos3Factor$BMC, Datos3Factor$PesoIntermedio, Datos3Factor$GPSurface,
              color = rainbow(length(levels(Datos3Factor$Dosis)))[Datos3Factor$Dosis], pch = 19,
              grid = TRUE, xlab = "BMC", ylab = "Peso Intermedio", zlab = "GPSurface",
              angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
       legend = levels(Datos3Factor$Dosis), fill = rainbow(length(levels(Datos3Factor$Dosis))))

Tratamiento de outliers

Todo el análisis que hemos hecho anteriormente ha sido sin quitar los outliers. Hay tres formas de gestionar los outliers: dejarlos, eliminarlos y cambiarlos por la media.

Estudio 1

Eliminar los outliers

# Función para eliminar outliers de columnas específicas
eliminar_outliers_columnas <- function(datos, columnas) {
  # Inicializar el dataframe sin outliers
  datos_sin_outliers <- datos
  
  for (columna in columnas) {
    # Calcular los cuartiles y el IQR de la columna
    Q1 <- quantile(datos[[columna]], 0.25)
    Q3 <- quantile(datos[[columna]], 0.75)
    IQR <- Q3 - Q1
    
    # Definir los límites inferior y superior
    limite_inferior <- Q1 - 1.5 * IQR
    limite_superior <- Q3 + 1.5 * IQR
    
    # Filtrar las filas que están dentro de los límites
    datos_sin_outliers <- datos_sin_outliers[datos_sin_outliers[[columna]] >= limite_inferior & datos_sin_outliers[[columna]] <= limite_superior, ]
  }
  
  return(datos_sin_outliers)
}

columnas_a_usar <- c(6:15)
Datos1_sin_outliers <- eliminar_outliers_columnas(Datos1, columnas_a_usar)
#print(Datos1_sin_outliers)

Comprobemos gráficamente la eliminación de los outliers.

par(mfrow = c(1, 2))
boxplot(Datos1[,6:15],main="Outliers antes de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))
boxplot(Datos1_sin_outliers[,6:15],main="Outliers después de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

Estandarizamos.

par(mfrow = c(1, 2))
sca<-scale(Datos1[,6:15])

boxplot(sca,main="Outliers antes de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

sca_sin_outliers<-scale(Datos1_sin_outliers[,6:15])

boxplot(sca_sin_outliers,main="Outliers después de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

Cambiar los outliers por la media

outlier<-function(data,na.rm=T){
  H<-1.5*IQR(data) 
  data[data<quantile(data,0.25,na.rm = T)-H]<-NA 
  data[data>quantile(data,0.75, na.rm = T)+H]<-NA 
  data[is.na(data)]<-mean(data, na.rm = T) 
  H<-1.5*IQR(data)
  
  if (TRUE %in% (data<quantile(data,0.25,na.rm = T)-H) | 
      TRUE %in% (data>quantile(data,0.75,na.rm = T)+H))
    outlier(data) 
  else
    return(data) 
}
Datos1_con_media<-Datos1

Datos1_con_media[[6]]<-outlier(Datos1_con_media[[6]]) 
Datos1_con_media[[7]]<-outlier(Datos1_con_media[[7]])  
Datos1_con_media[[8]]<-outlier(Datos1_con_media[[8]]) 
Datos1_con_media[[9]]<-outlier(Datos1_con_media[[9]]) 
Datos1_con_media[[10]]<-outlier(Datos1_con_media[[10]]) 
Datos1_con_media[[11]]<-outlier(Datos1_con_media[[11]]) 
Datos1_con_media[[12]]<-outlier(Datos1_con_media[[12]]) 
Datos1_con_media[[13]]<-outlier(Datos1_con_media[[13]]) 
Datos1_con_media[[14]]<-outlier(Datos1_con_media[[14]]) 
Datos1_con_media[[15]]<-outlier(Datos1_con_media[[15]]) 

Comprobemos gráficamente el cambio de los outliers por la media.

par(mfrow = c(1, 2))
boxplot(Datos1[,6:15],main="Outliers antes de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))
boxplot(Datos1_con_media[,6:15],main="Outliers después de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

Estandarizamos.

par(mfrow = c(1, 2))
sca<-scale(Datos1[,6:15])

boxplot(sca,main="Outliers antes de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

sca_sin_outliers1<-scale(Datos1_con_media[,6:15])

boxplot(sca_sin_outliers,main="Outliers después de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

Estudio 2

Eliminar los outliers

# Función para eliminar outliers de columnas específicas
eliminar_outliers_columnas <- function(datos, columnas) {
  # Inicializar el dataframe sin outliers
  datos_sin_outliers <- datos
  
  for (columna in columnas) {
    # Calcular los cuartiles y el IQR de la columna
    Q1 <- quantile(datos[[columna]], 0.25)
    Q3 <- quantile(datos[[columna]], 0.75)
    IQR <- Q3 - Q1
    
    # Definir los límites inferior y superior
    limite_inferior <- Q1 - 1.5 * IQR
    limite_superior <- Q3 + 1.5 * IQR
    
    # Filtrar las filas que están dentro de los límites
    datos_sin_outliers <- datos_sin_outliers[datos_sin_outliers[[columna]] >= limite_inferior & datos_sin_outliers[[columna]] <= limite_superior, ]
  }
  
  return(datos_sin_outliers)
}

columnas_a_usar <- c(6:15)
Datos2.2_sin_outliers <- eliminar_outliers_columnas(Datos2.2, columnas_a_usar)
#print(Datos2.2_sin_outliers)

Comprobemos gráficamente la eliminación de los outliers.

par(mfrow = c(1, 2))
boxplot(Datos2.2[,6:15],main="Outliers antes de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))
boxplot(Datos2.2_sin_outliers[,6:15],main="Outliers después de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

Estandarizamos.

par(mfrow = c(1, 2))
sca<-scale(Datos2.2[,6:15])

boxplot(sca,main="Outliers antes de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

sca_sin_outliers<-scale(Datos2.2_sin_outliers[,6:15])

boxplot(sca_sin_outliers,main="Outliers después de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

Cambiar los outliers por la media

outlier<-function(data,na.rm=T){
  H<-1.5*IQR(data) 
  data[data<quantile(data,0.25,na.rm = T)-H]<-NA 
  data[data>quantile(data,0.75, na.rm = T)+H]<-NA 
  data[is.na(data)]<-mean(data, na.rm = T) 
  H<-1.5*IQR(data)
  
  if (TRUE %in% (data<quantile(data,0.25,na.rm = T)-H) | 
      TRUE %in% (data>quantile(data,0.75,na.rm = T)+H))
    outlier(data) 
  else
    return(data) 
}

Datos2.2_con_media<-Datos2.2

#Llamamos a la función de outliers para cada variable en la que hemos identificado un outlier
#Sería suficiente fijarnos en las gráficas donde hay outliers y no llamar a la función para todas las variables


Datos2.2_con_media[[6]]<-outlier(Datos2.2_con_media[[6]]) 
Datos2.2_con_media[[7]]<-outlier(Datos2.2_con_media[[7]])  
Datos2.2_con_media[[8]]<-outlier(Datos2.2_con_media[[8]]) 
Datos2.2_con_media[[9]]<-outlier(Datos2.2_con_media[[9]]) 
Datos2.2_con_media[[10]]<-outlier(Datos2.2_con_media[[10]]) 
Datos2.2_con_media[[11]]<-outlier(Datos2.2_con_media[[11]]) 
Datos2.2_con_media[[12]]<-outlier(Datos2.2_con_media[[12]]) 
Datos2.2_con_media[[13]]<-outlier(Datos2.2_con_media[[13]]) 
Datos2.2_con_media[[14]]<-outlier(Datos2.2_con_media[[14]]) 
Datos2.2_con_media[[15]]<-outlier(Datos2.2_con_media[[15]])

Comprobemos gráficamente el cambio de los outliers por la media.

par(mfrow = c(1, 2))
boxplot(Datos2.2[,6:15],main="Outliers antes de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))
boxplot(Datos2.2_con_media[,6:15],main="Outliers después de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

Estandarizamos.

par(mfrow = c(1, 2))
sca<-scale(Datos2.2[,6:15])

boxplot(sca,main="Outliers antes de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

sca_sin_outliers2<-scale(Datos2.2_con_media[,6:15])

boxplot(sca_sin_outliers,main="Outliers después de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

Estudio 3

Eliminar los outliers

# Función para eliminar outliers de columnas específicas
eliminar_outliers_columnas <- function(datos, columnas) {
  # Inicializar el dataframe sin outliers
  datos_sin_outliers <- datos
  
  for (columna in columnas) {
    # Calcular los cuartiles y el IQR de la columna
    Q1 <- quantile(datos[[columna]], 0.25)
    Q3 <- quantile(datos[[columna]], 0.75)
    IQR <- Q3 - Q1
    
    # Definir los límites inferior y superior
    limite_inferior <- Q1 - 1.5 * IQR
    limite_superior <- Q3 + 1.5 * IQR
    
    # Filtrar las filas que están dentro de los límites
    datos_sin_outliers <- datos_sin_outliers[datos_sin_outliers[[columna]] >= limite_inferior & datos_sin_outliers[[columna]] <= limite_superior, ]
  }
  
  return(datos_sin_outliers)
}

columnas_a_usar <- c(6:15)
Datos3_sin_outliers <- eliminar_outliers_columnas(Datos3, columnas_a_usar)
#print(Datos3_sin_outliers)

Comprobemos gráficamente la eliminación de los outliers.

par(mfrow = c(1, 2))
boxplot(Datos3[,6:15],main="Outliers antes de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))
boxplot(Datos3_sin_outliers[,6:15],main="Outliers después de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

Estandarizamos.

par(mfrow = c(1, 2))
sca<-scale(Datos3[,6:15])

boxplot(sca,main="Outliers antes de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

sca_sin_outliers<-scale(Datos3_sin_outliers[,6:15])

boxplot(sca_sin_outliers,main="Outliers después de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

Cambiar los outliers por la media

outlier<-function(data,na.rm=T){
  H<-1.5*IQR(data) 
  data[data<quantile(data,0.25,na.rm = T)-H]<-NA 
  data[data>quantile(data,0.75, na.rm = T)+H]<-NA 
  data[is.na(data)]<-mean(data, na.rm = T) 
  H<-1.5*IQR(data)
  
  if (TRUE %in% (data<quantile(data,0.25,na.rm = T)-H) | 
      TRUE %in% (data>quantile(data,0.75,na.rm = T)+H))
    outlier(data) 
  else
    return(data) 
}

Datos3_con_media<-Datos3
#Llamamos a la función de outliers para cada variable en la que hemos identificado un outlier
#Sería suficiente fijarnos en las gráficas donde hay outliers y no llamar a la función para todas las variables

Datos3_con_media[[6]]<-outlier(Datos3_con_media[[6]]) 
Datos3_con_media[[7]]<-outlier(Datos3_con_media[[7]])  
Datos3_con_media[[8]]<-outlier(Datos3_con_media[[8]]) 
Datos3_con_media[[9]]<-outlier(Datos3_con_media[[9]]) 
Datos3_con_media[[10]]<-outlier(Datos3_con_media[[10]]) 
Datos3_con_media[[11]]<-outlier(Datos3_con_media[[11]]) 
Datos3_con_media[[12]]<-outlier(Datos3_con_media[[12]]) 
Datos3_con_media[[13]]<-outlier(Datos3_con_media[[13]]) 
Datos3_con_media[[14]]<-outlier(Datos3_con_media[[14]]) 
Datos3_con_media[[15]]<-outlier(Datos3_con_media[[15]])

Comprobemos gráficamente el cambio de los outliers por la media.

par(mfrow = c(1, 2))
boxplot(Datos3[,6:15],main="Outliers antes de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))
boxplot(Datos3_con_media[,6:15],main="Outliers después de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

Estandarizamos.

par(mfrow = c(1, 2))
sca<-scale(Datos3[,6:15])

boxplot(sca,main="Outliers antes de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

sca_sin_outliers3<-scale(Datos3_con_media[,6:15])

boxplot(sca_sin_outliers,main="Outliers después de quitarlos",
        xlab="All explanatory variables",
        ylab="Values",
        col=rainbow(10))

Análisis de Componentes Principales

Ahora, ya sin valores perdidos (los hemos eliminado), ni outliers (los hemos sustituido), hacemos ACP (Análisis de Componentes Principales)

Estudio de la correlación de las variables entre sí

Estudio 1

matriz_correlacion1<-cor(Datos1[6:15])
matriz_correlacion1
##                    PesoFinal       Tb.Sp         Tb.N         BMC         BMD
## PesoFinal       1.0000000000  0.04899019 -0.076027758  0.88857725  0.67063992
## Tb.Sp           0.0489901901  1.00000000 -0.966650769  0.03772223 -0.07477822
## Tb.N           -0.0760277581 -0.96665077  1.000000000 -0.07076835  0.06502617
## BMC             0.8885772453  0.03772223 -0.070768347  1.00000000  0.79200092
## BMD             0.6706399171 -0.07477822  0.065026167  0.79200092  1.00000000
## Longitud        0.7615998736  0.22750028 -0.274548124  0.82815136  0.58419435
## GPSurface       0.7479469259  0.16287707 -0.190656487  0.73317089  0.57053427
## GPVolume        0.3992803056  0.03884584 -0.008551150  0.20278668  0.16125363
## PesoInicial    -0.0002011703  0.00913795  0.001712557  0.06261252  0.01841571
## PesoIntermedio  0.8614173499  0.11924157 -0.148865043  0.87152256  0.69088510
##                   Longitud    GPSurface    GPVolume   PesoInicial
## PesoFinal       0.76159987  0.747946926  0.39928031 -0.0002011703
## Tb.Sp           0.22750028  0.162877069  0.03884584  0.0091379499
## Tb.N           -0.27454812 -0.190656487 -0.00855115  0.0017125570
## BMC             0.82815136  0.733170887  0.20278668  0.0626125219
## BMD             0.58419435  0.570534274  0.16125363  0.0184157082
## Longitud        1.00000000  0.660639694  0.09404165  0.2078779916
## GPSurface       0.66063969  1.000000000  0.60361811  0.0034486662
## GPVolume        0.09404165  0.603618114  1.00000000 -0.1495347667
## PesoInicial     0.20787799  0.003448666 -0.14953477  1.0000000000
## PesoIntermedio  0.84301621  0.589832006  0.01870350  0.0512494835
##                PesoIntermedio
## PesoFinal          0.86141735
## Tb.Sp              0.11924157
## Tb.N              -0.14886504
## BMC                0.87152256
## BMD                0.69088510
## Longitud           0.84301621
## GPSurface          0.58983201
## GPVolume           0.01870350
## PesoInicial        0.05124948
## PesoIntermedio     1.00000000
det(matriz_correlacion1)
## [1] 1.311426e-05

Estudio 2

matriz_correlacion2<-cor(Datos2.2[6:15])
matriz_correlacion2
##                 PesoFinal        Tb.Sp        Tb.N          BMC         BMD
## PesoFinal       1.0000000  0.185923507 -0.21026340  0.928349174  0.81306413
## Tb.Sp           0.1859235  1.000000000 -0.96382803 -0.006849368 -0.07397760
## Tb.N           -0.2102634 -0.963828027  1.00000000 -0.010983558  0.07778725
## BMC             0.9283492 -0.006849368 -0.01098356  1.000000000  0.89959506
## BMD             0.8130641 -0.073977602  0.07778725  0.899595064  1.00000000
## Longitud        0.8933615  0.081962798 -0.09378014  0.901471339  0.73605621
## GPSurface       0.7018824  0.058565382 -0.04145290  0.731245013  0.54036037
## GPVolume        0.1106489 -0.039730766  0.07318216  0.215404295  0.05951269
## PesoInicial     0.4002418 -0.028100641  0.05167504  0.442310869  0.37574638
## PesoIntermedio  0.9502286  0.248434899 -0.26003587  0.922008644  0.84160299
##                   Longitud   GPSurface    GPVolume PesoInicial PesoIntermedio
## PesoFinal       0.89336148  0.70188245  0.11064891  0.40024185      0.9502286
## Tb.Sp           0.08196280  0.05856538 -0.03973077 -0.02810064      0.2484349
## Tb.N           -0.09378014 -0.04145290  0.07318216  0.05167504     -0.2600359
## BMC             0.90147134  0.73124501  0.21540430  0.44231087      0.9220086
## BMD             0.73605621  0.54036037  0.05951269  0.37574638      0.8416030
## Longitud        1.00000000  0.79971969  0.38438797  0.41533525      0.8214429
## GPSurface       0.79971969  1.00000000  0.64932125  0.10495674      0.6367105
## GPVolume        0.38438797  0.64932125  1.00000000 -0.09734002      0.0551226
## PesoInicial     0.41533525  0.10495674 -0.09734002  1.00000000      0.3858588
## PesoIntermedio  0.82144287  0.63671048  0.05512260  0.38585884      1.0000000
det(matriz_correlacion2)
## [1] 4.090502e-07

Estudio 3

matriz_correlacion3<-cor(Datos3[6:15])
matriz_correlacion3
##                  PesoFinal       Tb.Sp       Tb.N         BMC        BMD
## PesoFinal       1.00000000 -0.38603705  0.4068459  0.95790423  0.8590239
## Tb.Sp          -0.38603705  1.00000000 -0.9723104 -0.50357610 -0.6405643
## Tb.N            0.40684592 -0.97231036  1.0000000  0.52230337  0.6476864
## BMC             0.95790423 -0.50357610  0.5223034  1.00000000  0.9322018
## BMD             0.85902394 -0.64056432  0.6476864  0.93220177  1.0000000
## Longitud        0.93550991 -0.35940363  0.3696785  0.94323641  0.8033721
## GPSurface       0.93119435 -0.32527309  0.3357800  0.90237031  0.8331614
## GPVolume        0.07754126  0.18163892 -0.1831888  0.05429013  0.0442623
## PesoInicial     0.13557640 -0.09253724  0.1228525  0.16003595  0.2106779
## PesoIntermedio  0.96651518 -0.35501921  0.3724031  0.94203518  0.8475858
##                   Longitud  GPSurface    GPVolume PesoInicial PesoIntermedio
## PesoFinal       0.93550991  0.9311944  0.07754126  0.13557640     0.96651518
## Tb.Sp          -0.35940363 -0.3252731  0.18163892 -0.09253724    -0.35501921
## Tb.N            0.36967854  0.3357800 -0.18318880  0.12285251     0.37240306
## BMC             0.94323641  0.9023703  0.05429013  0.16003595     0.94203518
## BMD             0.80337207  0.8331614  0.04426230  0.21067786     0.84758578
## Longitud        1.00000000  0.8771301  0.02949129  0.14211526     0.92169706
## GPSurface       0.87713006  1.0000000  0.35807064  0.19751925     0.88068747
## GPVolume        0.02949129  0.3580706  1.00000000  0.25670125    -0.01902111
## PesoInicial     0.14211526  0.1975193  0.25670125  1.00000000     0.05779478
## PesoIntermedio  0.92169706  0.8806875 -0.01902111  0.05779478     1.00000000
det(matriz_correlacion3)
## [1] 1.816053e-08

Test de Bartlett

Estudio 1

#dim(sca_sin_outliers1)
cortest.bartlett(cor(sca_sin_outliers1),n=56)
## $chisq
## [1] 248.3985
## 
## $p.value
## [1] 6.153321e-30
## 
## $df
## [1] 45

Estudio 2

#dim(sca_sin_outliers2)
cortest.bartlett(cor(sca_sin_outliers2),n=38)
## $chisq
## [1] 84.19442
## 
## $p.value
## [1] 0.0003579567
## 
## $df
## [1] 45

Estudio 3

#dim(sca_sin_outliers3)
cortest.bartlett(cor(sca_sin_outliers3),n=36)
## $chisq
## [1] 442.8086
## 
## $p.value
## [1] 8.5876e-67
## 
## $df
## [1] 45

Análisis de Componentes Principales

Estudio 1

PCA_1<-prcomp(Datos1_con_media[6:15], scale=T, center = T)
PCA_1$rotation
##                        PC1         PC2         PC3          PC4          PC5
## PesoFinal      -0.42237782 -0.01406992  0.18411363  0.306817447 -0.233732962
## Tb.Sp          -0.17948406 -0.56845073 -0.35663205 -0.002435807 -0.006438170
## Tb.N            0.11605409  0.62190190  0.22929194 -0.020449368 -0.169758706
## BMC            -0.43226109  0.19875483  0.01692992  0.111281614  0.009513443
## BMD            -0.33313043  0.14521134  0.06551204  0.168566068  0.875819166
## Longitud       -0.41128716 -0.02452259 -0.04627114  0.297477703 -0.273843241
## GPSurface      -0.36371262 -0.21247245  0.34235604 -0.256417443 -0.159664452
## GPVolume       -0.04739347 -0.25004961  0.65239142 -0.465139815  0.106184507
## PesoInicial    -0.20799929  0.15643373 -0.46055366 -0.632033255  0.072553849
## PesoIntermedio -0.36212499  0.31076004 -0.15712924 -0.307533504 -0.179826720
##                        PC6          PC7         PC8         PC9        PC10
## PesoFinal       0.19230758 -0.070989722 -0.69511730 -0.21604990  0.25854548
## Tb.Sp          -0.24792066 -0.243179428  0.12935726  0.21981729  0.57598751
## Tb.N           -0.46797906  0.001188083  0.05570919  0.12103088  0.53028359
## BMC            -0.04776174 -0.583724134  0.44868347 -0.44850427 -0.12408929
## BMD            -0.07092137  0.138155692 -0.07161163  0.18693243  0.06311439
## Longitud        0.01301219  0.695216366  0.42168371 -0.04309907  0.03804733
## GPSurface      -0.57986155 -0.014458407 -0.14486596  0.25849730 -0.43573252
## GPVolume        0.30788229  0.082702525  0.21554090 -0.16834823  0.32489087
## PesoInicial    -0.13438586  0.248168413 -0.20638466 -0.43514547  0.06190064
## PesoIntermedio  0.47537135 -0.155047874  0.05171300  0.60653394 -0.01068467

Vamos a calcular la desviación estándar de cada componente principal, la proporción de varianza explicada y la proporción de varianza explicada acumulada.

summary(PCA_1)
## Importance of components:
##                           PC1    PC2    PC3    PC4     PC5     PC6     PC7
## Standard deviation     1.9148 1.3710 1.1440 1.0408 0.80626 0.65430 0.61660
## Proportion of Variance 0.3666 0.1880 0.1309 0.1083 0.06501 0.04281 0.03802
## Cumulative Proportion  0.3666 0.5546 0.6855 0.7938 0.85881 0.90162 0.93964
##                            PC8     PC9    PC10
## Standard deviation     0.51334 0.45242 0.36792
## Proportion of Variance 0.02635 0.02047 0.01354
## Cumulative Proportion  0.96600 0.98646 1.00000

La siguiente gráfica muestra la proporción de varianza explicada.

Explained_variance <- PCA_1$sdev^2 / sum(PCA_1$sdev^2)

ggplot(data = data.frame(Explained_variance, pc = 1:10),
           aes(x = pc, y = Explained_variance, fill=Explained_variance )) + 
  geom_col(width = 0.3) + scale_y_continuous(limits = c(0,0.6)) + theme_bw() + 
  labs(x = "Principal component", y= "Proportion of variance")

La siguiente gráfica muestra la proporción de varianza explicada.

Cummulative_variance<-cumsum(Explained_variance)
ggplot( data = data.frame(Cummulative_variance, pc = 1:10),
            aes(x = pc, y = Cummulative_variance ,fill=Cummulative_variance )) + 
  geom_col(width = 0.5) + scale_y_continuous(limits = c(0,1)) + theme_bw() + 
  labs(x = "Principal component",
       y = "Proportion of cumulative variance")

Vemos cuántas componentes principales son adecuadas, usando la regla de Abdi: Las varianzas explicadas por los componentes principales se promedian, y se seleccionan aquellas cuya proporción de varianza explicada supera la media.

PCA_1$sdev^2
##  [1] 3.6664117 1.8796730 1.3087902 1.0832021 0.6500572 0.4281034 0.3801927
##  [8] 0.2635212 0.2046870 0.1353615
mean(PCA_1$sdev^2)
## [1] 1

Observamos que hay 4 por encima de la media, así que estudiamos el caso de 4 componentes principales. Esta gráfica muestra la proyección de las variables en dos dimensiones. Muestra el peso de la variable en la dirección de la componente principal.

# Gráfico 1: PC1 vs PC2
p1 <- fviz_pca_var(PCA_1, axes = c(1, 2), repel = TRUE, col.var = "cos2",
                   legend.title = "Distance", title = "PC1 vs PC2") + theme_bw()

# Gráfico 2: PC1 vs PC3
p2 <- fviz_pca_var(PCA_1, axes = c(1, 3), repel = TRUE, col.var = "cos2",
                   legend.title = "Distance", title = "PC1 vs PC3") + theme_bw()

# Gráfico 3: PC1 vs PC4
p3 <- fviz_pca_var(PCA_1, axes = c(1, 4), repel = TRUE, col.var = "cos2",
                   legend.title = "Distance", title = "PC1 vs PC4") + theme_bw()

# Gráfico 4: PC2 vs PC3
p4 <- fviz_pca_var(PCA_1, axes = c(2, 3), repel = TRUE, col.var = "cos2",
                   legend.title = "Distance", title = "PC2 vs PC3") + theme_bw()

# Gráfico 5: PC2 vs PC4
p5 <- fviz_pca_var(PCA_1, axes = c(2, 4), repel = TRUE, col.var = "cos2",
                   legend.title = "Distance", title = "PC2 vs PC4") + theme_bw()

# Gráfico 6: PC3 vs PC4
p6 <- fviz_pca_var(PCA_1, axes = c(3, 4), repel = TRUE, col.var = "cos2",
                   legend.title = "Distance", title = "PC3 vs PC4") + theme_bw()
grid.arrange(p1, p2, p3, p4, p5, p6, nrow = 3)

Estos resultados gráficos muestran las observaciones con su contribución a la varianza. Así como identificar con colores aquellas observaciones que explican la mayor varianza de los componentes principales.

# Gráfico 1: PC1 vs PC2
p1 <- fviz_pca_ind(PCA_1, col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
                   repel = TRUE, legend.title = "Contrib.var", title = "PC1 vs PC2") + theme_bw()

# Gráfico 2: PC1 vs PC3
p2 <- fviz_pca_ind(PCA_1, axes = c(1, 3), col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
                   repel = TRUE, legend.title = "Contrib.var", title = "PC1 vs PC3") + theme_bw()

# Gráfico 3: PC1 vs PC4
p3 <- fviz_pca_ind(PCA_1, axes = c(1, 4), col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
                   repel = TRUE, legend.title = "Contrib.var", title = "PC1 vs PC4") + theme_bw()

# Gráfico 4: PC2 vs PC3
p4 <- fviz_pca_ind(PCA_1, axes = c(2, 3), col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
                   repel = TRUE, legend.title = "Contrib.var", title = "PC2 vs PC3") + theme_bw()

# Gráfico 5: PC2 vs PC4
p5 <- fviz_pca_ind(PCA_1, axes = c(2, 4), col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
                   repel = TRUE, legend.title = "Contrib.var", title = "PC2 vs PC4") + theme_bw()

# Gráfico 6: PC3 vs PC4
p6 <- fviz_pca_ind(PCA_1, axes = c(3, 4), col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
                   repel = TRUE, legend.title = "Contrib.var", title = "PC3 vs PC4") + theme_bw()
grid.arrange(p1, p2, p3, p4, p5, p6, nrow = 3)

Estos resultados gráficos muestran las variables y las observaciones con su contribución a la varianza. Identifica las posibles relaciones entre las contribuciones de los registros a las varianzas de los componentes y el peso de las variables en cada componente principal.

# Gráfico 1: PC1 vs PC2
p1 <- fviz_pca(PCA_1, alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
               gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()

# Gráfico 2: PC1 vs PC3
p2 <- fviz_pca(PCA_1, axes = c(1, 3), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
               gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()

# Gráfico 3: PC1 vs PC4
p3 <- fviz_pca(PCA_1, axes = c(1, 4), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
               gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()

# Gráfico 4: PC2 vs PC3
p4 <- fviz_pca(PCA_1, axes = c(2, 3), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
               gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()

# Gráfico 5: PC2 vs PC4
p5 <- fviz_pca(PCA_1, axes = c(2, 4), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
               gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()

# Gráfico 6: PC3 vs PC4
p6 <- fviz_pca(PCA_1, axes = c(3, 4), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
               gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()

grid.arrange(p1, p2, p3, p4, p5, p6, nrow = 3)

Mostramos las coordenadas en el nuevo sistema de referencia

head(PCA_1$x)
##             PC1         PC2          PC3        PC4        PC5         PC6
## [1,]  0.2675163 -0.85141102 -0.155623192 -0.2840225 -1.2024258  0.27122398
## [2,]  2.8905478 -0.69771253 -0.117428533  0.9660797 -0.4560795 -0.84059860
## [3,] -2.3224893 -0.01339986  0.789150281 -0.6229904  0.5789967  0.29995198
## [4,] -1.3750855  1.37159909  0.008612358 -1.4599798 -1.1117157  0.93864317
## [5,] -3.0686746  0.89428864 -0.612825908 -0.5264444  0.1946402  0.04308362
## [6,]  1.4514801 -0.15087641  0.986361367 -0.1413722 -0.6230207 -0.30652335
##             PC7         PC8         PC9        PC10
## [1,] -0.1082483 -0.06499945 -0.80442331  0.48813868
## [2,] -0.5972964  0.20512951 -0.35830257 -0.26738714
## [3,] -1.3578986  0.01939389  0.02074547  0.27015198
## [4,]  0.1181536 -0.47237701  0.55526997  0.31625801
## [5,] -0.6622062  0.13964267 -0.59058134 -0.07965825
## [6,] -0.2994882  0.39283662  0.08906888 -0.10630711

Estudio 2

PCA_2<-prcomp(Datos2.2_con_media[6:15], scale=T, center = T)
PCA_2$rotation
##                        PC1        PC2         PC3         PC4         PC5
## PesoFinal      -0.25992492  0.2980351 -0.34757661 -0.47859323  0.17107702
## Tb.Sp          -0.23710824 -0.3098604  0.28815744 -0.40296686  0.43989353
## Tb.N            0.06134366  0.4273067 -0.42280589  0.35261251 -0.01130212
## BMC             0.33185093  0.1541500 -0.08552962  0.11718384  0.78389868
## BMD            -0.50447365 -0.1280594 -0.19453870  0.33384119 -0.01191529
## Longitud        0.33859316  0.2399386 -0.24789866 -0.53977919 -0.24429972
## GPSurface      -0.14483379 -0.4270193 -0.47112369 -0.22240881 -0.12003660
## GPVolume        0.12791088 -0.3892311 -0.52250284  0.11078784  0.24968223
## PesoInicial    -0.43493454  0.1285923 -0.10430309  0.01505287  0.04648235
## PesoIntermedio -0.40827143  0.4294838  0.07021918 -0.07563641  0.15448793
##                        PC6         PC7          PC8         PC9         PC10
## PesoFinal       0.22304724  0.02845598 -0.395331819  0.41670553  0.289819325
## Tb.Sp           0.08193385  0.44118628  0.451740845  0.04821205 -0.033018387
## Tb.N            0.14110777  0.44296051  0.483238167 -0.02806575  0.244919554
## BMC            -0.13025436 -0.44156481  0.124754826  0.05528321  0.001598339
## BMD             0.17456712 -0.25400188  0.213287489  0.46060872 -0.476473841
## Longitud       -0.17866716 -0.04959764  0.291311656  0.06432760 -0.540281481
## GPSurface       0.06163855 -0.41281942  0.304729359 -0.34523956  0.352953462
## GPVolume       -0.02092199  0.39980001 -0.400720968 -0.20776327 -0.347120891
## PesoInicial    -0.87666584  0.07513788 -0.009934938 -0.04509557  0.068619429
## PesoIntermedio  0.26736399 -0.10231783 -0.081670411 -0.66294812 -0.293627855

Vamos a calcular la desviación estándar de cada componente principal, la proporción de varianza explicada y la proporción de varianza explicada acumulada.

summary(PCA_2)
## Importance of components:
##                           PC1    PC2    PC3    PC4     PC5     PC6     PC7
## Standard deviation     1.5052 1.3747 1.2400 1.1146 0.96014 0.81996 0.68561
## Proportion of Variance 0.2266 0.1890 0.1538 0.1242 0.09219 0.06723 0.04701
## Cumulative Proportion  0.2266 0.4156 0.5693 0.6935 0.78574 0.85297 0.89998
##                            PC8     PC9    PC10
## Standard deviation     0.65022 0.61002 0.45311
## Proportion of Variance 0.04228 0.03721 0.02053
## Cumulative Proportion  0.94226 0.97947 1.00000

La siguiente gráfica muestra la proporción de varianza explicada.

Explained_variance <- PCA_2$sdev^2 / sum(PCA_2$sdev^2)

ggplot(data = data.frame(Explained_variance, pc = 1:10),
           aes(x = pc, y = Explained_variance, fill=Explained_variance )) + 
  geom_col(width = 0.3) + scale_y_continuous(limits = c(0,0.6)) + theme_bw() + 
  labs(x = "Principal component", y= "Proportion of variance")

La siguiente gráfica muestra la proporción de varianza explicada.

Cummulative_variance<-cumsum(Explained_variance)
ggplot( data = data.frame(Cummulative_variance, pc = 1:10),
            aes(x = pc, y = Cummulative_variance ,fill=Cummulative_variance )) + 
  geom_col(width = 0.5) + scale_y_continuous(limits = c(0,1)) + theme_bw() + 
  labs(x = "Principal component",
       y = "Proportion of cumulative variance")

Vemos cuántas componentes principales son adecuadas, usando la regla de Abdi: Las varianzas explicadas por los componentes principales se promedian, y se seleccionan aquellas cuya proporción de varianza explicada supera la media.

PCA_2$sdev^2
##  [1] 2.2657583 1.8898238 1.5376064 1.2423200 0.9218780 0.6723380 0.4700644
##  [8] 0.4227817 0.3721201 0.2053092
mean(PCA_2$sdev^2)
## [1] 1

Observamos que hay 4 por encima de la media, así que estudiamos el caso de 4 componentes principales. Esta gráfica muestra la proyección de las variables en dos dimensiones. Muestra el peso de la variable en la dirección de la componente principal.

# Gráfico 1: PC1 vs PC2
p1 <- fviz_pca_var(PCA_2, axes = c(1, 2), repel = TRUE, col.var = "cos2",
                   legend.title = "Distance", title = "PC1 vs PC2") + theme_bw()

# Gráfico 2: PC1 vs PC3
p2 <- fviz_pca_var(PCA_2, axes = c(1, 3), repel = TRUE, col.var = "cos2",
                   legend.title = "Distance", title = "PC1 vs PC3") + theme_bw()

# Gráfico 3: PC1 vs PC4
p3 <- fviz_pca_var(PCA_2, axes = c(1, 4), repel = TRUE, col.var = "cos2",
                   legend.title = "Distance", title = "PC1 vs PC4") + theme_bw()

# Gráfico 4: PC2 vs PC3
p4 <- fviz_pca_var(PCA_2, axes = c(2, 3), repel = TRUE, col.var = "cos2",
                   legend.title = "Distance", title = "PC2 vs PC3") + theme_bw()

# Gráfico 5: PC2 vs PC4
p5 <- fviz_pca_var(PCA_2, axes = c(2, 4), repel = TRUE, col.var = "cos2",
                   legend.title = "Distance", title = "PC2 vs PC4") + theme_bw()

# Gráfico 6: PC3 vs PC4
p6 <- fviz_pca_var(PCA_2, axes = c(3, 4), repel = TRUE, col.var = "cos2",
                   legend.title = "Distance", title = "PC3 vs PC4") + theme_bw()
grid.arrange(p1, p2, p3, p4, p5, p6, nrow = 3)

Estos resultados gráficos muestran las observaciones con su contribución a la varianza. Así como identificar con colores aquellas observaciones que explican la mayor varianza de los componentes principales.

# Gráfico 1: PC1 vs PC2
p1 <- fviz_pca_ind(PCA_2, col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
                   repel = TRUE, legend.title = "Contrib.var", title = "PC1 vs PC2") + theme_bw()

# Gráfico 2: PC1 vs PC3
p2 <- fviz_pca_ind(PCA_2, axes = c(1, 3), col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
                   repel = TRUE, legend.title = "Contrib.var", title = "PC1 vs PC3") + theme_bw()

# Gráfico 3: PC1 vs PC4
p3 <- fviz_pca_ind(PCA_2, axes = c(1, 4), col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
                   repel = TRUE, legend.title = "Contrib.var", title = "PC1 vs PC4") + theme_bw()

# Gráfico 4: PC2 vs PC3
p4 <- fviz_pca_ind(PCA_2, axes = c(2, 3), col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
                   repel = TRUE, legend.title = "Contrib.var", title = "PC2 vs PC3") + theme_bw()

# Gráfico 5: PC2 vs PC4
p5 <- fviz_pca_ind(PCA_2, axes = c(2, 4), col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
                   repel = TRUE, legend.title = "Contrib.var", title = "PC2 vs PC4") + theme_bw()

# Gráfico 6: PC3 vs PC4
p6 <- fviz_pca_ind(PCA_2, axes = c(3, 4), col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
                   repel = TRUE, legend.title = "Contrib.var", title = "PC3 vs PC4") + theme_bw()
grid.arrange(p1, p2, p3, p4, p5, p6, nrow = 3)

Estos resultados gráficos muestran las variables y las observaciones con su contribución a la varianza. Identifica las posibles relaciones entre las contribuciones de los registros a las varianzas de los componentes y el peso de las variables en cada componente principal.

# Gráfico 1: PC1 vs PC2
p1 <- fviz_pca(PCA_2, alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
               gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()

# Gráfico 2: PC1 vs PC3
p2 <- fviz_pca(PCA_2, axes = c(1, 3), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
               gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()

# Gráfico 3: PC1 vs PC4
p3 <- fviz_pca(PCA_2, axes = c(1, 4), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
               gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()

# Gráfico 4: PC2 vs PC3
p4 <- fviz_pca(PCA_2, axes = c(2, 3), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
               gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()

# Gráfico 5: PC2 vs PC4
p5 <- fviz_pca(PCA_2, axes = c(2, 4), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
               gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()

# Gráfico 6: PC3 vs PC4
p6 <- fviz_pca(PCA_2, axes = c(3, 4), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
               gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()

grid.arrange(p1, p2, p3, p4, p5, p6, nrow = 3)

Mostramos las coordenadas en el nuevo sistema de referencia

head(PCA_2$x)
##             PC1        PC2        PC3         PC4         PC5        PC6
## [1,] -0.5565609  1.0555072 -0.5221437 -0.14240038 -0.03484642  0.9084108
## [2,]  0.9333002 -0.9384008  2.7251929 -0.07547032 -0.89002717  0.4401783
## [3,]  2.2820594  0.6146465  2.1675107 -0.56366156 -0.12542127 -0.2913732
## [4,] -1.0898128  1.4132247 -1.9715946  0.53085606  0.96370299  0.6789528
## [5,] -0.5830088  1.2057056  1.0839525 -0.89863163  0.49375772  0.1674253
## [6,] -0.5097828 -0.9369064  0.5602292 -2.34375552 -1.54888116  0.4859421
##              PC7         PC8        PC9       PC10
## [1,] -0.46745970 -0.31353621 -0.5008955  0.0313710
## [2,] -0.30552079 -0.06875839 -1.1897007  0.9964603
## [3,] -0.63626742  0.34002247  0.1058022 -0.7520804
## [4,] -0.04304059  1.12403097 -0.8339104  0.6856212
## [5,]  0.33700562 -0.01534292 -0.1676559  0.2327255
## [6,]  1.56762171  0.76554558 -0.3195342 -0.3117803

Estudio 3

PCA_3<-prcomp(Datos3_con_media[6:15], scale=T, center = T)
PCA_3$rotation
##                        PC1         PC2         PC3         PC4          PC5
## PesoFinal      -0.39673358 -0.15725522 -0.05974109 -0.02209423  0.223708700
## Tb.Sp           0.25415886 -0.54303826 -0.22917115 -0.10898013  0.269494012
## Tb.N           -0.25208263  0.54499545  0.19743935  0.24929268 -0.003595032
## BMC            -0.40920080 -0.06996250 -0.02757952  0.02344307  0.195991225
## BMD            -0.39659089  0.05588512  0.04864212  0.12142489 -0.039844026
## Longitud       -0.38263075 -0.16813439 -0.12425103 -0.07929067  0.365364278
## GPSurface      -0.37830852 -0.26824077  0.17112292  0.01737851  0.046277032
## GPVolume       -0.02017475 -0.42080586  0.79776660 -0.02411437 -0.296473049
## PesoInicial     0.14501040 -0.21612759 -0.03968514  0.94702119  0.123818798
## PesoIntermedio -0.28066312 -0.22656417 -0.46805425  0.07858178 -0.773598333
##                         PC6          PC7         PC8          PC9         PC10
## PesoFinal       0.228149851 -0.235384793 -0.49970271 -0.522602781 -0.365370643
## Tb.Sp          -0.439670214  0.437303981 -0.33921780  0.008999114  0.061605549
## Tb.N           -0.067634099  0.624420509 -0.37609358  0.007990931  0.048027822
## BMC            -0.147175330  0.001963112  0.28414183 -0.418533383  0.714399034
## BMD            -0.776490179 -0.253618142  0.13139464  0.143887971 -0.340099181
## Longitud        0.284816556  0.417927267  0.49451461  0.221286421 -0.347580355
## GPSurface       0.175315256 -0.213181226 -0.36457661  0.664063133  0.321131282
## GPVolume        0.004249656  0.192082291  0.11074756 -0.196659068 -0.098755359
## PesoInicial     0.090358748 -0.059573312  0.07651613 -0.016323749 -0.025900613
## PesoIntermedio  0.073669048  0.196413387 -0.01415723 -0.044483427  0.002842317

Vamos a calcular la desviación estándar de cada componente principal, la proporción de varianza explicada y la proporción de varianza explicada acumulada.

summary(PCA_3)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5    PC6     PC7
## Standard deviation     2.3934 1.2892 1.0217 0.94026 0.61297 0.3377 0.31351
## Proportion of Variance 0.5728 0.1662 0.1044 0.08841 0.03757 0.0114 0.00983
## Cumulative Proportion  0.5728 0.7390 0.8434 0.93181 0.96938 0.9808 0.99061
##                            PC8     PC9    PC10
## Standard deviation     0.23525 0.16921 0.09942
## Proportion of Variance 0.00553 0.00286 0.00099
## Cumulative Proportion  0.99615 0.99901 1.00000

La siguiente gráfica muestra la proporción de varianza explicada.

Explained_variance <- PCA_3$sdev^2 / sum(PCA_3$sdev^2)

ggplot(data = data.frame(Explained_variance, pc = 1:10),
           aes(x = pc, y = Explained_variance, fill=Explained_variance )) + 
  geom_col(width = 0.3) + scale_y_continuous(limits = c(0,0.6)) + theme_bw() + 
  labs(x = "Principal component", y= "Proportion of variance")

La siguiente gráfica muestra la proporción de varianza explicada.

Cummulative_variance<-cumsum(Explained_variance)
ggplot( data = data.frame(Cummulative_variance, pc = 1:10),
            aes(x = pc, y = Cummulative_variance ,fill=Cummulative_variance )) + 
  geom_col(width = 0.5) + scale_y_continuous(limits = c(0,1)) + theme_bw() + 
  labs(x = "Principal component",
       y = "Proportion of cumulative variance")

Vemos cuántas componentes principales son adecuadas, usando la regla de Abdi: Las varianzas explicadas por los componentes principales se promedian, y se seleccionan aquellas cuya proporción de varianza explicada supera la media.

PCA_3$sdev^2
##  [1] 5.728162672 1.662065418 1.043793886 0.884082735 0.375734284 0.114015776
##  [7] 0.098286304 0.055342812 0.028632712 0.009883402
mean(PCA_3$sdev^2)
## [1] 1

Observamos que hay 3 por encima de la media, así que estudiamos el caso de 3 componentes principales. Esta gráfica muestra la proyección de las variables en dos dimensiones. Muestra el peso de la variable en la dirección de la componente principal.

# Gráfico 1: PC1 vs PC2
p1 <- fviz_pca_var(PCA_3, axes = c(1, 2), repel = TRUE, col.var = "cos2",
                   legend.title = "Distance", title = "PC1 vs PC2") + theme_bw()

# Gráfico 2: PC1 vs PC3
p2 <- fviz_pca_var(PCA_3, axes = c(1, 3), repel = TRUE, col.var = "cos2",
                   legend.title = "Distance", title = "PC1 vs PC3") + theme_bw()

# Gráfico 3: PC2 vs PC3
p3 <- fviz_pca_var(PCA_3, axes = c(2, 3), repel = TRUE, col.var = "cos2",
                   legend.title = "Distance", title = "PC2 vs PC3") + theme_bw()

grid.arrange(p1, p2, p3, nrow = 2)

Estos resultados gráficos muestran las observaciones con su contribución a la varianza. Así como identificar con colores aquellas observaciones que explican la mayor varianza de los componentes principales.

# Gráfico 1: PC1 vs PC2
p1<-fviz_pca_ind(PCA_3,col.ind = "contrib",
                 gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), 
                 repel=TRUE,legend.title="Contrib.var", title="Records")+theme_bw()

# Gráfico 2: PC1 vs PC3
p2<-fviz_pca_ind(PCA_3,axes=c(1,3),col.ind = "contrib", 
                 gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
                 repel=TRUE,legend.title="Contrib.var", title="Records")+theme_bw()

# Gráfico 1: PC2 vs PC3
p3<-fviz_pca_ind(PCA_3,axes=c(2,3),col.ind = "contrib", 
                 gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
                 repel=TRUE,legend.title="Contrib.var", title="Records")+theme_bw()

grid.arrange(p1, p2, p3, nrow = 2)

Estos resultados gráficos muestran las variables y las observaciones con su contribución a la varianza. Identifica las posibles relaciones entre las contribuciones de los registros a las varianzas de los componentes y el peso de las variables en cada componente principal.

# Gráfico 1: PC1 vs PC2
p1 <- fviz_pca(PCA_3, alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
               gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()

# Gráfico 2: PC1 vs PC3
p2 <- fviz_pca(PCA_3, axes = c(1, 3), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
               gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()

# Gráfico 3: PC2 vs PC3
p3 <- fviz_pca(PCA_3, axes = c(2, 3), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
               gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()

grid.arrange(p1, p2, p3, nrow = 2)

Mostramos las coordenadas en el nuevo sistema de referencia

head(PCA_3$x)
##             PC1        PC2        PC3         PC4         PC5         PC6
## [1,] -1.5233677 -0.3581962  0.9210322  0.41030042  0.12089330  0.02626068
## [2,]  3.1185773  0.6346817 -0.5727801  0.03208267 -0.96206638 -0.33121675
## [3,]  1.0295563  0.4364019  1.4661299  0.54437498  0.01583101  0.29865500
## [4,]  0.5514543  0.9841047  0.3842833 -0.81213671 -0.66415803 -0.25097263
## [5,]  1.2337493  0.9170182  0.3097037  0.75243837  0.16870026  0.39776375
## [6,]  0.4911084  1.0929372 -0.4006754  0.16191980 -0.08825771  0.25773892
##              PC7         PC8         PC9        PC10
## [1,]  0.06277742  0.09463181  0.18848194 -0.15272820
## [2,]  0.20938983  0.22941818 -0.25554403 -0.07244920
## [3,] -0.22116901  0.02723566  0.13598678  0.05547987
## [4,]  0.16641356  0.12817490  0.02373183  0.02365107
## [5,] -0.05732605  0.02661690  0.26384796 -0.09626547
## [6,] -0.01416743 -0.21702662 -0.02839555  0.05045322
#print(PCA_3$x)
#print(PCA_3)

Estudio de la normalidad

Previo a la construcción de métodos de clasificación se debe analizar la normalidad multivariante de los datos con el test propuesto en clase de prácticas de análisis discriminante.

Estudio 1

PCA_1scores<-as.data.frame(PCA_1$x[,1:4])
outliers1 <- mvn(data = PCA_1scores, mvnTest = "hz", multivariateOutlierMethod = "quan")

Estudio 2

PCA_2scores<-as.data.frame(PCA_2$x[,1:4])
outliers2 <- mvn(data = PCA_2scores, mvnTest = "hz", multivariateOutlierMethod = "quan")

Estudio 3

PCA_3scores<-as.data.frame(PCA_3$x[,1:3])
outliers3 <- mvn(data = PCA_3scores, mvnTest = "hz", multivariateOutlierMethod = "quan")

Royston multivariate normality test

Estudio 1

royston_test1 <- mvn(data = PCA_1scores, mvnTest = "royston", multivariatePlot = "qq")

royston_test1$multivariateNormality
##      Test        H     p value MVN
## 1 Royston 17.61511 0.001467199  NO

Estudio 2

royston_test2 <- mvn(data = PCA_2scores, mvnTest = "royston", multivariatePlot = "qq")

royston_test2$multivariateNormality
##      Test        H   p value MVN
## 1 Royston 2.532671 0.6387953 YES

Estudio 3

royston_test3 <- mvn(data = PCA_3scores, mvnTest = "royston", multivariatePlot = "qq")

royston_test3$multivariateNormality
##      Test        H    p value MVN
## 1 Royston 10.02739 0.01833478  NO

Henze-Zirkler multivariate normality test

Estudio 1

hz_test1 <- mvn(data = PCA_1scores, mvnTest = "hz") 
hz_test1$multivariateNormality
##            Test       HZ    p value MVN
## 1 Henze-Zirkler 0.929499 0.07517825 YES

Estudio 2

hz_test2 <- mvn(data = PCA_2scores, mvnTest = "hz") 
hz_test2$multivariateNormality
##            Test       HZ   p value MVN
## 1 Henze-Zirkler 0.811583 0.2443798 YES

Estudio 3

hz_test3 <- mvn(data = PCA_3scores, mvnTest = "hz") 
hz_test3$multivariateNormality
##            Test       HZ     p value MVN
## 1 Henze-Zirkler 1.246543 0.000942336  NO

Las variables respuesta (dependientes) serían: Tb.Sp, Tb.N, BMC, BMD, Longitud, GPSurface, GPVolume.

Las variables independientes (predictoras) serían: Grupo, Dosis, DiasDosis,Peso, PesoInicial, Peso Intermedio, DiasEspera.

Análisis clúster

Reescalamos las componentes principales.

PCA_1scores_scale<-PCA_1scores
PCA_1scores_scale<-as.data.frame(scale(PCA_1scores_scale))

PCA_2scores_scale<-PCA_2scores
PCA_2scores_scale<-as.data.frame(scale(PCA_2scores_scale))

PCA_3scores_scale<-PCA_3scores
PCA_3scores_scale<-as.data.frame(scale(PCA_3scores_scale))

Hacemos la matriz de distancia para cada estudio.

Estudio 1

distance<- get_dist(PCA_1scores_scale)
fviz_dist(distance, gradient = list(low ="yellow", mid = "white", high = "brown"))

Estudio 2

distance<- get_dist(PCA_2scores_scale)
fviz_dist(distance, gradient = list(low ="yellow", mid = "white", high = "brown"))

Estudio 3

distance<- get_dist(PCA_3scores_scale)
fviz_dist(distance, gradient = list(low ="yellow", mid = "white", high = "brown"))

Clustering jerárquico - método de Ward

Hacemos los dendogramas.

Estudio 1

dendrogram <- hclust(dist(PCA_1scores_scale, method = 'euclidean'), method = 'ward.D') 
ggdendrogram(dendrogram, rotate = FALSE, labels = FALSE, theme_dendro = TRUE) +
  labs(title = "Dendrograma")

Estudio 2

dendrogram <- hclust(dist(PCA_2scores_scale, method = 'euclidean'), method = 'ward.D') 
ggdendrogram(dendrogram, rotate = FALSE, labels = FALSE, theme_dendro = TRUE) +
  labs(title = "Dendrograma")

Estudio 3

dendrogram <- hclust(dist(PCA_3scores_scale, method = 'euclidean'), method = 'ward.D') 
ggdendrogram(dendrogram, rotate = FALSE, labels = FALSE, theme_dendro = TRUE) +
  labs(title = "Dendrograma")

Clustering no jerárquico - algoritmo K-means

Estudio 1

set.seed(123)

k2_1 <- kmeans(PCA_1scores_scale, centers = 2, nstart = 25)
k3_1 <- kmeans(PCA_1scores_scale, centers = 3, nstart = 25) 
k4_1 <- kmeans(PCA_1scores_scale, centers = 4, nstart = 25) 
k5_1 <- kmeans(PCA_1scores_scale, centers = 5, nstart = 25)
k6_1 <- kmeans(PCA_1scores_scale, centers = 6, nstart = 25)
k7_1 <- kmeans(PCA_1scores_scale, centers = 7, nstart = 25)
k9_1 <- kmeans(PCA_1scores_scale, centers = 9, nstart = 25)
# Plots to compare
fviz_cluster(k2_1, data = PCA_1scores_scale) + ggtitle("k = 2") 

fviz_cluster(k3_1, data = PCA_1scores_scale) + ggtitle("k = 3") 

fviz_cluster(k4_1, data = PCA_1scores_scale) + ggtitle("k = 4") 

fviz_cluster(k5_1, data = PCA_1scores_scale) + ggtitle("k = 5")

fviz_cluster(k6_1, data = PCA_1scores_scale) + ggtitle("k = 6") 

fviz_cluster(k7_1, data = PCA_1scores_scale) + ggtitle("k = 7")

fviz_cluster(k9_1, data = PCA_1scores_scale) + ggtitle("k = 9") 

Apliquemos el método de Elbow, de Silhouette y el del estadístico de brecha para ver el número óptimo de clústers.

Método de Elbow
set.seed(123)
fviz_nbclust(PCA_1scores_scale, kmeans, method = "wss")

Método de Silhouette
set.seed(123)
fviz_nbclust(PCA_1scores_scale, kmeans, method = "silhouette")

Método estadístico de brecha (GAP)
set.seed(123)
gap_stat <- clusGap(PCA_1scores_scale, FUN = kmeans, nstart = 25,K.max = 10, B = 50) 
fviz_gap_stat(gap_stat)

Estudio 2

set.seed(123)
k2_2 <- kmeans(PCA_2scores_scale, centers = 2, nstart = 25)
k3_2 <- kmeans(PCA_2scores_scale, centers = 3, nstart = 25) 
k4_2 <- kmeans(PCA_2scores_scale, centers = 4, nstart = 25) 
k5_2 <- kmeans(PCA_2scores_scale, centers = 5, nstart = 25)
# Plots to compare
fviz_cluster(k2_2, data = PCA_2scores_scale) + ggtitle("k = 2") 

fviz_cluster(k3_2, data = PCA_2scores_scale) + ggtitle("k = 3") 

fviz_cluster(k4_2, data = PCA_2scores_scale) + ggtitle("k = 4") 

fviz_cluster(k5_2, data = PCA_2scores_scale) + ggtitle("k = 5")

Apliquemos el método de Elbow, de Silhouette y el del estadístico de brecha para ver el número óptimo de clústers.

Método de Elbow
set.seed(123)
fviz_nbclust(PCA_2scores_scale, kmeans, method = "wss")

Método de Silhouette
set.seed(123)
fviz_nbclust(PCA_2scores_scale, kmeans, method = "silhouette")

Método estadístico de brecha (GAP)
set.seed(123)
gap_stat <- clusGap(PCA_2scores_scale, FUN = kmeans, nstart = 25,K.max = 10, B = 50) 
fviz_gap_stat(gap_stat)

Estudio 3 ( sin escalar )

set.seed(123)
k2_3 <- kmeans(PCA_3scores, centers = 2, nstart = 25)
k3_3 <- kmeans(PCA_3scores, centers = 3, nstart = 25) 
k4_3 <- kmeans(PCA_3scores, centers = 4, nstart = 25) 
k5_3 <- kmeans(PCA_3scores, centers = 5, nstart = 25)
# Plots to compare
fviz_cluster(k2_3, data = PCA_3scores) + ggtitle("k = 2")

fviz_cluster(k3_3, data = PCA_3scores) + ggtitle("k = 3")

fviz_cluster(k4_3, data = PCA_3scores) + ggtitle("k = 4")

fviz_cluster(k5_3, data = PCA_3scores) + ggtitle("k = 5")

Apliquemos el método de Elbow, de Silhouette y el del estadístico de brecha para ver el número óptimo de clústers.

Método de Elbow
set.seed(123)
fviz_nbclust(PCA_3scores_scale, kmeans, method = "wss")

Método de Silhouette
set.seed(123)
fviz_nbclust(PCA_3scores_scale, kmeans, method = "silhouette")

Método estadístico de brecha (GAP)
set.seed(123)
gap_stat <- clusGap(PCA_3scores_scale, FUN = kmeans, nstart = 25,K.max = 10, B = 50) 
fviz_gap_stat(gap_stat)

Cálculo de vectores de media

Veamos el número óptimo de clusters. Calculamos vectores de medias de 2, 3 y 4 clusters. #### Estudio 1 ##### Vectores de medias con 2 cluster.

Datos1_con_media$Cluster2 <- k2_1$cluster
Datos1_con_media$Cluster3 <- k3_1$cluster
Datos1_con_media$Cluster4 <- k4_1$cluster

colMeans(subset(Datos1_con_media[3:17], Cluster2 == 1))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      1.0000000     21.1000000     14.0000000    288.5552830      0.3245529 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      2.9881437      0.2632654    142.2262667     38.7577212     39.1877675 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      4.7023645    123.1900000    164.6585354     14.0000000      1.0000000
colMeans(subset(Datos1_con_media[3:17], Cluster2 == 2))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      1.0000000     36.8846154     14.0000000    304.2852685      0.2933026 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      3.2657038      0.2682001    143.4288177     38.9815385     41.2414331 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.9588533    116.5961538    165.1677350     14.0000000      2.0000000
Vectores de medias con 3 cluster.
colMeans(subset(Datos1_con_media[3:18], Cluster3 == 1))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      1.0000000     25.7600000     14.0000000    277.7183396      0.3170226 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      2.9973644      0.2509190    138.9594800     38.2376655     38.5942069 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.0396320    118.5800000    158.3110480     14.0000000      1.2000000 
##       Cluster3 
##      1.0000000
colMeans(subset(Datos1_con_media[3:18], Cluster3 == 2))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      1.0000000     36.4615385     14.0000000    303.2198839      0.2967179 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      3.3144923      0.2766495    145.1298718     39.0961538     41.5172458 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.7030430    122.4230769    171.1128739     14.0000000      1.7307692 
##       Cluster3 
##      2.0000000
colMeans(subset(Datos1_con_media[3:18], Cluster3 == 3))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      1.0000000      0.0000000     14.0000000    348.2800000      0.3444445 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      2.6883400      0.2810602    149.7147188     40.7620000     40.7213439 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      4.3462400    115.9400000    165.4812500     14.0000000      1.4000000 
##       Cluster3 
##      3.0000000
colMeans(subset(Datos1_con_media[3:18], Dosis == 80))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      1.0000000     80.0000000     14.0000000    285.6750000      0.2910465 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      3.3513250      0.2647625    142.1288750     38.7700000     41.3159225 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.6107819    119.7375000    162.1545139     14.0000000      1.6250000 
##       Cluster3 
##      1.7500000
Vectores de medias con 4 cluster.
colMeans(subset(Datos1_con_media[3:19], Cluster4 == 1))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      1.0000000     21.7000000     14.0000000    275.2250000      0.3052372 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      3.1449455      0.2545940    139.7934000     38.3180000     37.8528443 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      4.6114900    120.5550000    161.5147475     14.0000000      1.0000000 
##       Cluster3       Cluster4 
##      1.1000000      1.0000000
colMeans(subset(Datos1_con_media[3:19], Cluster4 == 2))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      1.0000000     39.4210526     14.0000000    295.2587885      0.2854012 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      3.4086947      0.2624125    141.9950877     38.6578947     40.8136171 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      6.0205033    115.6631579    161.9559211     14.0000000      2.0000000 
##       Cluster3       Cluster4 
##      1.7368421      2.0000000
colMeans(subset(Datos1_con_media[3:19], Cluster4 == 3))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      1.0000000      0.0000000     14.0000000    346.6250000      0.3304307 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      2.7367000      0.2849361    149.0038985     40.7675000     40.0302841 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      4.3432750    109.7250000    165.4812500     14.0000000      1.5000000 
##       Cluster3       Cluster4 
##      3.0000000      3.0000000
colMeans(subset(Datos1_con_media[3:19], Cluster4 == 4))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##       1.000000      31.461538      14.000000     312.858345       0.347182 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##       2.764746       0.281054     146.626692      39.409357      42.713349 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##       5.539127     129.200000     174.210363      14.000000       1.384615 
##       Cluster3       Cluster4 
##       1.923077       4.000000

Estudio 2

Vectores de medias con 2 cluster.
Datos2.2_con_media$Cluster2 <- k2_2$cluster
Datos2.2_con_media$Cluster3 <- k3_2$cluster
Datos2.2_con_media$Cluster4 <- k4_2$cluster

colMeans(subset(Datos2.2_con_media[3:17], Cluster2 == 1))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      2.0000000    176.8421053     21.0000000    351.7894737      0.3148076 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      3.1184895      0.3245605    154.3637895     41.7489522     43.2144266 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.4325263    115.9684211    164.3157895     28.0000000      1.0000000
colMeans(subset(Datos2.2_con_media[3:17], Cluster2 == 2))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      2.0000000     58.9473684     21.0000000    331.9638596      0.3190193 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      2.8985947      0.3242006    164.3222632     41.7358057     43.8724149 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.8739789    116.2157895    154.9175439     28.0000000      2.0000000
Vectores de medias con 3 cluster.
colMeans(subset(Datos2.2_con_media[3:18], Cluster3 == 1))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      2.0000000    175.0000000     21.0000000    356.4312500      0.3167109 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      3.1983813      0.3244072    156.0120000     41.7480842     43.2287417 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.4088875    117.3312500    169.1937500     28.0000000      1.0000000 
##       Cluster3 
##      1.0000000
colMeans(subset(Datos2.2_con_media[3:18], Cluster3 == 2))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      2.0000000      8.8888889     21.0000000    341.0237037      0.3253313 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      2.7106444      0.3234977    176.6955556     41.7344344     45.6307988 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.9622333    120.6555556    161.0037037     28.0000000      2.0000000 
##       Cluster3 
##      2.0000000
colMeans(subset(Datos2.2_con_media[3:18], Cluster3 == 3))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##       2.000000     123.076923      21.000000     324.553846       0.311335 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##       2.981131       0.324959     151.429462      41.740857      42.485610 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##       5.740100     111.407692     146.869231      28.000000       1.769231 
##       Cluster3 
##       3.000000
Vectores de medias con 4 cluster.
colMeans(subset(Datos2.2_con_media[3:19], Cluster4 == 1))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      2.0000000      8.8888889     21.0000000    341.0237037      0.3253313 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      2.7106444      0.3234977    176.6955556     41.7344344     45.6307988 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.9622333    120.6555556    161.0037037     28.0000000      2.0000000 
##       Cluster3       Cluster4 
##      2.0000000      1.0000000
colMeans(subset(Datos2.2_con_media[3:19], Cluster4 == 2))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      2.0000000    160.0000000     21.0000000    354.6500000      0.3189563 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      3.1786250      0.3242581    157.4050000     41.7466503     42.7150195 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.0384500    117.8583333    175.7666667     28.0000000      1.0000000 
##       Cluster3       Cluster4 
##      1.0000000      2.0000000
colMeans(subset(Datos2.2_con_media[3:19], Cluster4 == 3))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      2.0000000     70.0000000     21.0000000    321.5750000      0.3133542 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      2.8665125      0.3247514    151.8255000     41.7371055     41.3465789 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.1177875    111.9750000    151.6625000     28.0000000      1.8750000 
##       Cluster3       Cluster4 
##      3.0000000      3.0000000
colMeans(subset(Datos2.2_con_media[3:19], Cluster4 == 4))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      2.0000000    213.3333333     21.0000000    343.7444444      0.3089356 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      3.2059111      0.3250970    151.2567778     41.7493159     44.5133260 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      6.6399778    112.8333333    143.7666667     28.0000000      1.3333333 
##       Cluster3       Cluster4 
##      2.1111111      4.0000000
colMeans(subset(Datos2.2_con_media[3:19], Dosis == 80))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      2.0000000     80.0000000     21.0000000    339.7090909      0.3234777 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      3.1073455      0.3247939    153.7318182     41.7446471     44.1898289 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.8196273    113.6636364    158.6727273     28.0000000      1.3636364 
##       Cluster3       Cluster4 
##      2.0000000      2.8181818

Estudio 3

Vectores de medias con 2 cluster.
Datos3_con_media$Cluster2 <- k2_3$cluster
Datos3_con_media$Cluster3 <- k3_3$cluster
Datos3_con_media$Cluster4 <- k4_3$cluster
Datos3_con_media$Cluster5 <- k5_3$cluster

colMeans(subset(Datos3_con_media[3:17], Cluster2 == 1))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      3.0000000    519.2307692     14.0000000    211.9230769      0.3209199 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      3.1546231      0.1793077    118.0496923     33.9794231     33.9432231 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.6273692     71.7447653     99.1657484     14.0000000      1.0000000
colMeans(subset(Datos3_con_media[3:17], Cluster2 == 2))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##        3.00000        0.00000       14.00000      306.80000        0.27303 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##        3.63439        0.26381      137.91720       37.09100       41.00512 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##        5.69381       71.71914       99.70932       14.00000        2.00000
Vectores de medias con 3 cluster.
colMeans(subset(Datos3_con_media[3:18], Cluster3 == 1))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##        3.00000        0.00000       14.00000      306.80000        0.27303 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##        3.63439        0.26381      137.91720       37.09100       41.00512 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##        5.69381       71.71914       99.70932       14.00000        2.00000 
##       Cluster3 
##        1.00000
colMeans(subset(Datos3_con_media[3:18], Cluster3 == 2))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      3.0000000    230.7692308     14.0000000    216.6307692      0.2685615 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      3.6757538      0.1907077    122.8738462     34.1626923     34.4485692 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.5732538     71.7269413     99.1175006     14.0000000      1.0000000 
##       Cluster3 
##      2.0000000
colMeans(subset(Datos3_con_media[3:18], Cluster3 == 3))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      3.0000000    807.6923077     14.0000000    207.2153846      0.3732782 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      2.6334923      0.1679077    113.2255385     33.7961538     33.4378769 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.6814846     71.7625893     99.2139962     14.0000000      1.0000000 
##       Cluster3 
##      3.0000000
Vectores de medias con 4 cluster.
colMeans(subset(Datos3_con_media[3:19], Cluster4 == 1))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##       3.000000       0.000000      14.000000     312.183333       0.251250 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##       3.868433       0.268300     138.779333      37.171667      40.009517 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##       4.855683      71.701187      99.725402      14.000000       2.000000 
##       Cluster3       Cluster4 
##       1.000000       1.000000
colMeans(subset(Datos3_con_media[3:19], Cluster4 == 2))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      3.0000000    230.7692308     14.0000000    216.6307692      0.2685615 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      3.6757538      0.1907077    122.8738462     34.1626923     34.4485692 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.5732538     71.7269413     99.1175006     14.0000000      1.0000000 
##       Cluster3       Cluster4 
##      2.0000000      2.0000000
colMeans(subset(Datos3_con_media[3:19], Cluster4 == 3))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      3.0000000    807.6923077     14.0000000    207.2153846      0.3732782 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      2.6334923      0.1679077    113.2255385     33.7961538     33.4378769 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.6814846     71.7625893     99.2139962     14.0000000      1.0000000 
##       Cluster3       Cluster4 
##      3.0000000      3.0000000
colMeans(subset(Datos3_con_media[3:19], Cluster4 == 4))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##       3.000000       0.000000      14.000000     298.725000       0.305700 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##       3.283325       0.257075     136.624000      36.970000      42.498525 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##       6.951000      71.746067      99.685195      14.000000       2.000000 
##       Cluster3       Cluster4 
##       1.000000       4.000000

Conclusión número de clusters

Estudio 1

Decidimos elegir k=3. Por tanto, los datos se quedan agrupados en 3 clusters.

set.seed(123)
k3_1 <- kmeans(PCA_1scores_scale, centers = 3, nstart = 25) 
fviz_cluster(k3_1, data = PCA_1scores_scale) + ggtitle("k = 3") 

Vector de medias de cada cluster:

colMeans(subset(Datos1_con_media[3:18], Cluster3 == 1))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      1.0000000     25.7600000     14.0000000    277.7183396      0.3170226 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      2.9973644      0.2509190    138.9594800     38.2376655     38.5942069 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.0396320    118.5800000    158.3110480     14.0000000      1.2000000 
##       Cluster3 
##      1.0000000
colMeans(subset(Datos1_con_media[3:18], Cluster3 == 2))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      1.0000000     36.4615385     14.0000000    303.2198839      0.2967179 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      3.3144923      0.2766495    145.1298718     39.0961538     41.5172458 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.7030430    122.4230769    171.1128739     14.0000000      1.7307692 
##       Cluster3 
##      2.0000000
colMeans(subset(Datos1_con_media[3:18], Cluster3 == 3))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      1.0000000      0.0000000     14.0000000    348.2800000      0.3444445 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      2.6883400      0.2810602    149.7147188     40.7620000     40.7213439 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      4.3462400    115.9400000    165.4812500     14.0000000      1.4000000 
##       Cluster3 
##      3.0000000

Estudio 2

Decidimos elegir k=4. Por tanto, los datos se quedan agrupados en 4 clusters.

set.seed(123)
k4_2 <- kmeans(PCA_2scores_scale, centers = 4, nstart = 25) 
fviz_cluster(k4_2, data = PCA_2scores_scale) + ggtitle("k = 4") 

Vector de medias de cada cluster:

colMeans(subset(Datos2.2_con_media[3:19], Cluster4 == 1))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      2.0000000      8.8888889     21.0000000    341.0237037      0.3253313 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      2.7106444      0.3234977    176.6955556     41.7344344     45.6307988 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.9622333    120.6555556    161.0037037     28.0000000      2.0000000 
##       Cluster3       Cluster4 
##      2.0000000      1.0000000
colMeans(subset(Datos2.2_con_media[3:19], Cluster4 == 2))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      2.0000000    160.0000000     21.0000000    354.6500000      0.3189563 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      3.1786250      0.3242581    157.4050000     41.7466503     42.7150195 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.0384500    117.8583333    175.7666667     28.0000000      1.0000000 
##       Cluster3       Cluster4 
##      1.0000000      2.0000000
colMeans(subset(Datos2.2_con_media[3:19], Cluster4 == 3))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      2.0000000     70.0000000     21.0000000    321.5750000      0.3133542 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      2.8665125      0.3247514    151.8255000     41.7371055     41.3465789 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.1177875    111.9750000    151.6625000     28.0000000      1.8750000 
##       Cluster3       Cluster4 
##      3.0000000      3.0000000
colMeans(subset(Datos2.2_con_media[3:19], Cluster4 == 4))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      2.0000000    213.3333333     21.0000000    343.7444444      0.3089356 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      3.2059111      0.3250970    151.2567778     41.7493159     44.5133260 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      6.6399778    112.8333333    143.7666667     28.0000000      1.3333333 
##       Cluster3       Cluster4 
##      2.1111111      4.0000000

Estudio 3

Decidimos elegir k=3. Por tanto, los datos se quedan agrupados en 3 clusters.

set.seed(123)
k3_3 <- kmeans(PCA_3scores, centers = 3, nstart = 25)
fviz_cluster(k3_3, data = PCA_3scores) + ggtitle("k = 3")

Vector de medias de cada cluster:

colMeans(subset(Datos3_con_media[3:18], Cluster3 == 1))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##        3.00000        0.00000       14.00000      306.80000        0.27303 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##        3.63439        0.26381      137.91720       37.09100       41.00512 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##        5.69381       71.71914       99.70932       14.00000        2.00000 
##       Cluster3 
##        1.00000
colMeans(subset(Datos3_con_media[3:18], Cluster3 == 2))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      3.0000000    230.7692308     14.0000000    216.6307692      0.2685615 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      3.6757538      0.1907077    122.8738462     34.1626923     34.4485692 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.5732538     71.7269413     99.1175006     14.0000000      1.0000000 
##       Cluster3 
##      2.0000000
colMeans(subset(Datos3_con_media[3:18], Cluster3 == 3))
##        Estudio          Dosis      DiasDosis      PesoFinal          Tb.Sp 
##      3.0000000    807.6923077     14.0000000    207.2153846      0.3732782 
##           Tb.N            BMC            BMD       Longitud      GPSurface 
##      2.6334923      0.1679077    113.2255385     33.7961538     33.4378769 
##       GPVolume    PesoInicial PesoIntermedio     DiasEspera       Cluster2 
##      5.6814846     71.7625893     99.2139962     14.0000000      1.0000000 
##       Cluster3 
##      3.0000000

ANÁLISIS DISCRIMINANTE LINEAL

ESTUDIO 1

set.seed(123)
PCA_1scores_scale$Cluster3 <- k3_1$cluster
datos<-PCA_1scores_scale[,-4]
datos$Cluster3 <- as.factor(datos$Cluster3)
set.seed(123)
trainIndex<-createDataPartition(datos$Cluster3,p=0.80)$Resample1 
datos_train<-datos[trainIndex,] 
datos_test<-datos[-trainIndex,]

Nuestras variables predictoras van a ser las 3 primeras componentes principales y la dosis.

plot(datos_train[, 1:3], col = datos$Cluster3, pch = 19, main="Datos de entrenamiento")

plot(datos_test[, 1:3], col = datos$Cluster3, pch = 19, main="Datos de test")

En primer lugar exploramos cómo de bien (o mal) clasifica en los grupos de efectividad de la dosis cada una de las variables explicativas consideradas de forma independiente. Para ello dibujamos los histogramas superpuestos. Si los histogramas se separan, la variable considerada sería un buen clasificador individual para el grupo.

p1 <- ggplot(data = datos_train, aes(x = PC1, fill = Cluster3)) + geom_histogram(position = "identity", alpha = 0.5)
p2 <- ggplot(data = datos_train, aes(x = PC2, fill = Cluster3)) + geom_histogram(position = "identity", alpha = 0.5)
p3 <- ggplot(data = datos_train, aes(x = PC3, fill = Cluster3)) + geom_histogram(position = "identity", alpha = 0.5)
ggarrange(p1, p2, p3, nrow = 3, common.legend = TRUE)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Parece que la variable PC3 de la tercera componente principal es la que mejor separa entre los dos grupos.

Veamos qué pares de variables separan mejor:

pairs(x = datos_train[, c("PC1","PC2","PC3")], 
      col = c("blue", "orange")[datos_train$Cluster3], pch = 19)

Las combinaciones PC2-PC3 y PC1-PC3 parecen separar bien los datos.

Veamos como de bien separan las tres variables y si por tanto tiene sentido usarlas para hacer el análisis discriminante.

scatterplot3d(datos_train$PC1, datos_train$PC2, datos_train$PC3, 
              color = c("blue", "orange")[datos_train$Cluster3], pch = 19,
              grid = TRUE, xlab = "var_exp1", ylab = "var_exp2",
              zlab = "var_exp3", angle = 65, cex.axis = 0.6) 
legend("topleft",
         bty = "n", cex = 0.8,
         title = "efecividad_dosis",
         c("a", "b"), fill = c("blue", "orange"))

Separan bastante bien, así que tiene sentido proceder con el análisis discriminante.

A continuación hacemos una exploración gráfica de la normalidad de las distribuciones univariantes de nuestros predictores representando los histogramas y los gráficos qqplots.

par(mfcol = c(2, 3)) 
for (k in 1:3) {
  j0 <- names(datos_train)[k]
  x0 <- seq(min(datos_train[, k]), max(datos_train[, k]), le = 50) 
  for (i in 1:2) {
    i0 <- levels(datos_train$Cluster3)[i]
    x <- datos[datos_train$Cluster3 == i0, j0]
    hist(x, proba = T, col = grey(0.8), main = paste("efectividad_dosis", i0), xlab = j0) 
    lines(x0, dnorm(x0, mean(x), sd(x)), col = "blue", lwd = 2)
  } 
}

##### Representation of normal quantiles of each variable for each species
par(mfrow=c(2,3)) 
for (k in 1:3) {
  j0 <- names(datos_train)[k]
  x0 <- seq(min(datos_train[, k]), max(datos_train[, k]), le = 50) 
  for (i in 1:2) {
    i0 <- levels(datos_train$Cluster3)[i]
    x <- datos[datos_train$Cluster3 == i0, j0]
    qqnorm(x, main = paste("efectividad_dosis", i0, j0), pch = 19, col = i + 1) 
    qqline(x)
  } 
}

#####Normalidad univariante
datos_tidy <- melt(datos_train, value.name = "valor") 
## Using Cluster3 as id variables
datos_tidy %>%
  group_by(Cluster3,variable) %>%
  summarise(p_value_Shapiro.test = round(shapiro.test(valor)$p.value,5))
## `summarise()` has grouped output by 'Cluster3'. You can override using the
## `.groups` argument.
## # A tibble: 9 × 3
## # Groups:   Cluster3 [3]
##   Cluster3 variable p_value_Shapiro.test
##   <fct>    <fct>                   <dbl>
## 1 1        PC1                    0.590 
## 2 1        PC2                    0.0597
## 3 1        PC3                    0.136 
## 4 2        PC1                    0.700 
## 5 2        PC2                    0.370 
## 6 2        PC3                    0.438 
## 7 3        PC1                    0.243 
## 8 3        PC2                    0.524 
## 9 3        PC3                    0.284

Observamos que el p-valor>0.05, podemos asumir normalidad univariante en cada variable.

Tests de normalidad multivariante:

outliers <- mvn(data = datos_train[,-4], mvnTest = "hz", multivariateOutlierMethod = "quan")

royston_test <- mvn(data = datos_train[,-4], mvnTest = "royston", multivariatePlot = "qq")

royston_test$multivariateNormality
##      Test        H    p value MVN
## 1 Royston 8.019127 0.04562044  NO
hz_test <- mvn(data = datos_train[,-4], mvnTest = "hz") 
hz_test$multivariateNormality
##            Test        HZ   p value MVN
## 1 Henze-Zirkler 0.6837849 0.3736335 YES

Homogeneidad de varianzas

boxM(data = datos_train[,-4], grouping = datos_train[, 4])
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  datos_train[, -4]
## Chi-Sq (approx.) = 13.061, df = 12, p-value = 0.3646

Observamos que se da normalidad y homogeneidad de las varianzas.

modelo_lda <- lda(formula = Cluster3 ~ PC1 + PC2 + PC3, data=datos_train)
modelo_lda
## Call:
## lda(Cluster3 ~ PC1 + PC2 + PC3, data = datos_train)
## 
## Prior probabilities of groups:
##          1          2          3 
## 0.08888889 0.46666667 0.44444444 
## 
## Group means:
##          PC1        PC2        PC3
## 1 -1.3718383 -0.4563999 -0.3148362
## 2 -0.5783738  0.3208546  0.4149847
## 3  0.8157435 -0.1246100 -0.4954859
## 
## Coefficients of linear discriminants:
##            LD1        LD2
## PC1  1.6060635 -0.4559855
## PC2 -0.6357729 -0.6652696
## PC3 -0.7826050 -0.6374750
## 
## Proportion of trace:
##    LD1    LD2 
## 0.9513 0.0487
#probabilidades: 0.6198, 0.3902

Una vez construido el clasificador, podemos clasificar nuevos datos en función de sus medidas sin más que llamar a la función predict.

Veamos para todos los datos test.

nuevas_observaciones<-datos_test
predict(object = modelo_lda, newdata = nuevas_observaciones) 
## $class
##  [1] 2 3 2 2 2 2 2 3 3 3 3
## Levels: 1 2 3
## 
## $posterior
##               1            2            3
## 4  6.372456e-02 9.347727e-01 0.0015027191
## 7  1.640030e-02 1.226638e-01 0.8609359134
## 13 1.605429e-02 9.835519e-01 0.0003938489
## 21 1.960784e-02 9.606043e-01 0.0197878602
## 24 2.714848e-02 9.721102e-01 0.0007413120
## 25 2.748075e-01 7.146342e-01 0.0105582776
## 33 4.689169e-02 7.266001e-01 0.2265082433
## 37 1.281882e-05 8.759697e-05 0.9998995842
## 42 3.187139e-04 4.470158e-03 0.9952111282
## 50 1.124855e-05 7.747008e-04 0.9992140507
## 52 6.855588e-04 1.749723e-03 0.9975647187
## 
## $x
##           LD1        LD2
## 4  -1.7567557 -0.3552411
## 7   0.7677199  0.5384874
## 13 -2.1031712 -1.4725951
## 21 -0.9112295 -1.1209042
## 24 -1.9357527 -1.0470107
## 25 -1.1737072  1.0225313
## 33 -0.1440404 -0.1452081
## 37  3.0223824  0.9563182
## 42  1.8586148  0.2403761
## 50  2.4826781 -0.8464948
## 52  2.0525918  1.5382039

Validación del modelo

pred <- predict(object = modelo_lda, newdata = datos_test) 
confusionmatrix(datos_test$Cluster3, pred$class)
##   new 1 new 2 new 3
## 1     0     1     0
## 2     0     5     0
## 3     0     0     5

Clasificación del porcentaje del error

trainig_error <- mean(datos_test$Cluster3 != pred$class) * 100 
paste("trainig_error=", trainig_error, "%")
## [1] "trainig_error= 9.09090909090909 %"

Visualización de las clasificaciones

partimat(Cluster3 ~ PC1 + PC2 + PC3, 
         data = datos_test, method = "lda", prec = 200,
         image.colors = c("green", "orange","blue"),
         col.mean = "yellow",nplots.vert =1, nplots.hor=3)

Training error: 9%. Cierto error en PC1-PC2, 18%; PC1-PC3, 9%.

ESTUDIO 2

PCA_2scores_scale$Cluster4 <- k4_2$cluster
datos<-PCA_2scores_scale[,-4]
datos$Cluster4 <- as.factor(datos$Cluster4)
# Partitioning the dataset: training (80%) + test (20%)
set.seed(123)
trainIndex<-createDataPartition(datos$Cluster4,p=0.80)$Resample1 
datos_train<-datos[trainIndex,] 
datos_test<-datos[-trainIndex,]

Nuestras variables predictoras van a ser las 3 primeras componentes principales y la dosis.

plot(datos_train[, 1:3], col = datos$Cluster3, pch = 19, main="Datos de entrenamiento")

plot(datos_test[, 1:3], col = datos$Cluster3, pch = 19, main="Datos de test")

Exploramos clasificación.

p1 <- ggplot(data = datos_train, aes(x = PC1, fill = Cluster4)) + geom_histogram(position = "identity", alpha = 0.5)
p2 <- ggplot(data = datos_train, aes(x = PC2, fill = Cluster4)) + geom_histogram(position = "identity", alpha = 0.5)
p3 <- ggplot(data = datos_train, aes(x = PC3, fill = Cluster4)) + geom_histogram(position = "identity", alpha = 0.5)
ggarrange(p1, p2, p3, nrow = 3, common.legend = TRUE)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Parece que la variable PC3 de la tercera componente principal es la que mejor separa entre los dos grupos.

Veamos qué pares de variables separan mejor:

pairs(x = datos_train[, c("PC1","PC2","PC3")], 
      col = c("blue", "orange")[datos_train$Cluster4], pch = 19)

Las combinaciones PC2-PC3 y PC1-PC3 parecen separar bien los datos

Veamos como de bien separan las tres variables y si por tanto tiene sentido usarlas para hacer el análisis discriminante

scatterplot3d(datos_train$PC1, datos_train$PC2, datos_train$PC3, 
              color = c("blue", "orange")[datos_train$Cluster4], pch = 19,
              grid = TRUE, xlab = "var_exp1", ylab = "var_exp2",
              zlab = "var_exp3", angle = 65, cex.axis = 0.6) 
legend("topleft",
       bty = "n", cex = 0.8,
       title = "efecividad_dosis",
       c("a", "b"), fill = c("blue", "orange"))

Separan bastante bien, así que tiene sentido proceder con el análisis discriminante.

A continuación hacemos una exploración gráfica de la normalidad de las distribuciones univariantes de nuestros predictores representando los histogramas y los gráficos qqplots.

par(mfcol = c(2, 3)) 
for (k in 1:3) {
  j0 <- names(datos_train)[k]
  x0 <- seq(min(datos_train[, k]), max(datos_train[, k]), le = 50) 
  for (i in 1:2) {
    i0 <- levels(datos_train$Cluster4)[i]
    x <- datos[datos_train$Cluster4 == i0, j0]
    hist(x, proba = T, col = grey(0.8), main = paste("efectividad_dosis", i0), xlab = j0) 
    lines(x0, dnorm(x0, mean(x), sd(x)), col = "blue", lwd = 2)
  } 
}

par(mfrow=c(2,3)) 
for (k in 1:3) {
  j0 <- names(datos_train)[k]
  x0 <- seq(min(datos_train[, k]), max(datos_train[, k]), le = 50) 
  for (i in 1:2) {
    i0 <- levels(datos_train$Cluster4)[i]
    x <- datos[datos_train$Cluster4 == i0, j0]
    qqnorm(x, main = paste("efectividad_dosis", i0, j0), pch = 19, col = i + 1) 
    qqline(x)
  } 
}

#Normalidad univariante
datos_tidy <- melt(datos_train, value.name = "valor") 
## Using Cluster4 as id variables
datos_tidy %>%
  group_by(Cluster4,variable) %>%
  summarise(p_value_Shapiro.test = round(shapiro.test(valor)$p.value,5))
## `summarise()` has grouped output by 'Cluster4'. You can override using the
## `.groups` argument.
## # A tibble: 12 × 3
## # Groups:   Cluster4 [4]
##    Cluster4 variable p_value_Shapiro.test
##    <fct>    <fct>                   <dbl>
##  1 1        PC1                     0.349
##  2 1        PC2                     0.641
##  3 1        PC3                     0.586
##  4 2        PC1                     0.146
##  5 2        PC2                     0.149
##  6 2        PC3                     0.766
##  7 3        PC1                     0.530
##  8 3        PC2                     0.562
##  9 3        PC3                     0.795
## 10 4        PC1                     0.813
## 11 4        PC2                     0.792
## 12 4        PC3                     0.367

Observamos que el p-valor>0.05, podemos asumir normalidad univariante en cada.

Tests de normalidad multivariante:

outliers <- mvn(data = datos_train[,-4], mvnTest = "hz", multivariateOutlierMethod = "quan")

royston_test <- mvn(data = datos_train[,-4], mvnTest = "royston", multivariatePlot = "qq")

royston_test$multivariateNormality
##      Test        H   p value MVN
## 1 Royston 1.998275 0.5727658 YES
hz_test <- mvn(data = datos_train[,-4], mvnTest = "hz") 
hz_test$multivariateNormality
##            Test        HZ   p value MVN
## 1 Henze-Zirkler 0.6472712 0.3817672 YES

Homogeneidad de varianzas

boxM(data = datos_train[,-4], grouping = datos_train[, 4])
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  datos_train[, -4]
## Chi-Sq (approx.) = 20.802, df = 18, p-value = 0.2895

Observamos que se da normalidad y homogeneidad de las varianzas.

modelo_lda <- lda(formula = Cluster4 ~ PC1 + PC2 + PC3, data=datos_train)
modelo_lda
## Call:
## lda(Cluster4 ~ PC1 + PC2 + PC3, data = datos_train)
## 
## Prior probabilities of groups:
##       1       2       3       4 
## 0.21875 0.25000 0.25000 0.28125 
## 
## Group means:
##          PC1         PC2        PC3
## 1  0.6599792  0.03286496  1.0466227
## 2 -0.5456697  0.91151269  0.4214908
## 3 -1.1064918 -1.02213585 -0.2658523
## 4  0.5758493  0.03657006 -1.0758341
## 
## Coefficients of linear discriminants:
##            LD1         LD2        LD3
## PC1 -1.3886224  0.97816263 -0.5633177
## PC2 -1.0444490 -0.04111225  0.9622804
## PC3 -0.8087304 -1.34473390 -0.4414437
## 
## Proportion of trace:
##    LD1    LD2    LD3 
## 0.5642 0.3511 0.0847
#probabilidades: 0.9513, 0.0487

Una vez construido el clasificador, podemos clasificar nuevos datos en función de sus medidas sin más que llamar a la función predict

Veamos para todos los datos test

nuevas_observaciones<-datos_test
predict(object = modelo_lda, newdata = nuevas_observaciones) 
## $class
## [1] 1 4 4 3 2 2
## Levels: 1 2 3 4
## 
## $posterior
##               1            2            3            4
## 2  0.9827994437 0.0171845072 5.716587e-06 1.033252e-05
## 11 0.3769804388 0.0025946913 8.099923e-05 6.203439e-01
## 19 0.0055702677 0.0003012557 2.306756e-07 9.941282e-01
## 21 0.0003815275 0.1069681004 8.917523e-01 8.980249e-04
## 33 0.3384897107 0.4389712756 7.762480e-07 2.225382e-01
## 38 0.2714866078 0.4313807576 4.669169e-03 2.924635e-01
## 
## $x
##           LD1        LD2         LD3
## 2  -2.1123001 -2.2635332 -2.04196261
## 11 -1.3713408  1.6426389 -2.26599200
## 19 -1.9791146  3.0688773 -0.39825960
## 21  1.8261149 -1.4142720  0.75168656
## 33 -2.7972879  0.9137039  1.49106874
## 38 -0.6681194  0.2290076 -0.03096111

Validación del modelo

pred <- predict(object = modelo_lda, newdata = datos_test) 
confusionmatrix(datos_test$Cluster4, pred$class)
##   new 1 new 2 new 3 new 4
## 1     1     0     0     0
## 2     0     2     0     0
## 3     0     0     1     0
## 4     0     0     0     2

Clasificación del porcentaje del error

trainig_error <- mean(datos_test$Cluster4 != pred$class) * 100 
paste("trainig_error=", trainig_error, "%")
## [1] "trainig_error= 0 %"

Visualización de las clasificaciones

partimat(Cluster4 ~ PC1 + PC2 + PC3, 
         data = datos_test, method = "lda", prec = 200,
         image.colors = c("green", "orange","blue","red"),
         col.mean = "yellow",nplots.vert =1, nplots.hor=3)

Training error: 0%. Rrror de 40% en PC1-PC3, sustancial por falta de datos.

ESTUDIO 3

PCA_3scores$Cluster3 <- k3_3$cluster
datos<-PCA_3scores
datos$Cluster3 <- as.factor(datos$Cluster3)
# Partitioning the dataset: training (80%) + test (20%)
set.seed(123)
trainIndex<-createDataPartition(datos$Cluster3,p=0.80)$Resample1 
datos_train<-datos[trainIndex,] 
datos_test<-datos[-trainIndex,]

Nuestras variables predictoras van a ser las 3 primeras componentes principales y la dosis.

plot(datos_train[, 1:3], col = datos$Cluster3, pch = 19, main="Datos de entrenamiento")

plot(datos_test[, 1:3], col = datos$Cluster3, pch = 19, main="Datos de test")

Exploramos clasificación.

p1 <- ggplot(data = datos_train, aes(x = PC1, fill = Cluster3)) + geom_histogram(position = "identity", alpha = 0.5)
p2 <- ggplot(data = datos_train, aes(x = PC2, fill = Cluster3)) + geom_histogram(position = "identity", alpha = 0.5)
p3 <- ggplot(data = datos_train, aes(x = PC3, fill = Cluster3)) + geom_histogram(position = "identity", alpha = 0.5)
ggarrange(p1, p2, p3, nrow = 3, common.legend = TRUE)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Parece que la variable PC3 de la tercera componente principal es la que mejor separa entre los dos grupos.

Veamos qué pares de variables separan mejor:

pairs(x = datos_train[, c("PC1","PC2","PC3")], 
      col = c("blue", "orange")[datos_train$Cluster3], pch = 19)

Las combinaciones PC2-PC3 y PC1-PC3 parecen separar bien los datos

Veamos como de bien separan las tres variables y si por tanto tiene sentido usarlas para hacer el análisis discriminante

scatterplot3d(datos_train$PC1, datos_train$PC2, datos_train$PC3, 
              color = c("blue", "orange")[datos_train$Cluster3], pch = 19,
              grid = TRUE, xlab = "var_exp1", ylab = "var_exp2",
              zlab = "var_exp3", angle = 65, cex.axis = 0.6) 
legend("topleft",
       bty = "n", cex = 0.8,
       title = "efecividad_dosis",
       c("a", "b"), fill = c("blue", "orange"))

Separan bastante bien, así que tiene sentido proceder con el análisis discriminante.

A continuación hacemos una exploración gráfica de la normalidad de las distribuciones univariantes de nuestros predictores representando los histogramas y los gráficos qqplots.

par(mfcol = c(2, 3)) 
for (k in 1:3) {
  j0 <- names(datos_train)[k]
  x0 <- seq(min(datos_train[, k]), max(datos_train[, k]), le = 50) 
  for (i in 1:2) {
    i0 <- levels(datos_train$Cluster3)[i]
    x <- datos[datos_train$Cluster3 == i0, j0]
    hist(x, proba = T, col = grey(0.8), main = paste("efectividad_dosis", i0), xlab = j0) 
    lines(x0, dnorm(x0, mean(x), sd(x)), col = "blue", lwd = 2)
  } 
}

par(mfrow=c(2,3)) 
for (k in 1:3) {
  j0 <- names(datos_train)[k]
  x0 <- seq(min(datos_train[, k]), max(datos_train[, k]), le = 50) 
  for (i in 1:2) {
    i0 <- levels(datos_train$Cluster3)[i]
    x <- datos[datos_train$Cluster3 == i0, j0]
    qqnorm(x, main = paste("efectividad_dosis", i0, j0), pch = 19, col = i + 1) 
    qqline(x)
  } 
}

#Normalidad univariante
datos_tidy <- melt(datos_train, value.name = "valor") 
## Using Cluster3 as id variables
datos_tidy %>%
  group_by(Cluster3,variable) %>%
  summarise(p_value_Shapiro.test = round(shapiro.test(valor)$p.value,5))
## `summarise()` has grouped output by 'Cluster3'. You can override using the
## `.groups` argument.
## # A tibble: 9 × 3
## # Groups:   Cluster3 [3]
##   Cluster3 variable p_value_Shapiro.test
##   <fct>    <fct>                   <dbl>
## 1 1        PC1                    0.110 
## 2 1        PC2                    0.348 
## 3 1        PC3                    0.0920
## 4 2        PC1                    0.0484
## 5 2        PC2                    0.0903
## 6 2        PC3                    0.257 
## 7 3        PC1                    0.217 
## 8 3        PC2                    0.799 
## 9 3        PC3                    0.212

Observamos que el p-valor>0.05, podemos asumir normalidad univariante en todas excepto en una.

Tests de normalidad multivariante:

outliers <- mvn(data = datos_train[,-4], mvnTest = "hz", multivariateOutlierMethod = "quan")

royston_test <- mvn(data = datos_train[,-4], mvnTest = "royston", multivariatePlot = "qq")

royston_test$multivariateNormality
##      Test        H    p value MVN
## 1 Royston 9.794117 0.02039988  NO
hz_test <- mvn(data = datos_train[,-4], mvnTest = "hz") 
hz_test$multivariateNormality
##            Test       HZ     p value MVN
## 1 Henze-Zirkler 1.175856 0.001754329  NO

Homogeneidad de varianzas

boxM(data = datos_train[,-4], grouping = datos_train[, 4])
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  datos_train[, -4]
## Chi-Sq (approx.) = 16.557, df = 12, p-value = 0.167

Observamos que se da normalidad y homogeneidad de las varianzas.

modelo_lda <- lda(formula = Cluster3 ~ PC1 + PC2 + PC3, data=datos_train)
modelo_lda
## Call:
## lda(Cluster3 ~ PC1 + PC2 + PC3, data = datos_train)
## 
## Prior probabilities of groups:
##         1         2         3 
## 0.3666667 0.2666667 0.3666667 
## 
## Group means:
##          PC1        PC2        PC3
## 1  2.0553136 -0.6432026 -0.3127310
## 2 -3.4443121 -0.4688611 -0.1386789
## 3  0.3644586  1.0726328  0.6388531
## 
## Coefficients of linear discriminants:
##           LD1         LD2
## PC1 1.2431245  0.05057212
## PC2 0.1505819 -1.07116237
## PC3 0.1192883 -0.72051310
## 
## Proportion of trace:
##    LD1    LD2 
## 0.8462 0.1538
#probabilidades:

Una vez construido el clasificador, podemos clasificar nuevos datos en función de sus medidas sin más que llamar a la función predict

Veamos para todos los datos test

nuevas_observaciones<-datos_test
predict(object = modelo_lda, newdata = nuevas_observaciones) 
## $class
## [1] 3 3 2 2 1 1
## Levels: 1 2 3
## 
## $posterior
##               1            2            3
## 4  2.139487e-02 1.761181e-07 9.786049e-01
## 7  1.632248e-02 7.016742e-11 9.836775e-01
## 10 4.180382e-10 9.999998e-01 1.807804e-07
## 16 2.524699e-10 9.999995e-01 4.887665e-07
## 30 9.999656e-01 1.401340e-15 3.441981e-05
## 32 9.997673e-01 3.234498e-11 2.326927e-04
## 
## $x
##           LD1        LD2
## 4   0.9036463 -1.2073013
## 7   2.0945594 -2.1026849
## 10 -4.2218932  1.3303718
## 16 -4.2477080  0.7730017
## 30  3.9206753  2.1962734
## 32  2.4179979  2.4601340

Validación del modelo

pred <- predict(object = modelo_lda, newdata = datos_test) 
confusionmatrix(datos_test$Cluster3, pred$class)
##   new 1 new 2 new 3
## 1     2     0     0
## 2     0     2     0
## 3     0     0     2

Clasificación del porcentaje del error

trainig_error <- mean(datos_test$Cluster3 != pred$class) * 100 
paste("trainig_error=", trainig_error, "%")
## [1] "trainig_error= 0 %"

Visualización de las clasificaciones

partimat(Cluster3 ~ PC1 + PC2 + PC3, 
         data = datos_test, method = "lda", prec = 200,
         image.colors = c("green", "orange","blue"),
         col.mean = "yellow",nplots.vert =1, nplots.hor=3)

Training error: 0%. Error en PC2-PC3 del 33’3%.